Pan-Cancer Analysis for Identification of Tumor Antigens and Immune Subtypes in mRNA Vaccine Development

Methods

We collected normalized gene expression data and clinical data from the Chinese Glioma Genome Atlas (CGGA) and The Cancer Genome Atlas (TCGA). Gene Expression Profiling Interactive Analysis (GEPIA) version 2 was used to analyze the differential expression of cancer genes and to calculate the prognostic index of each selected antigen. We used the cBio Cancer Genomics Portal (cBioPortal) to visualize the gene alterations of potential tumor antigens in the TCGA data and the Tumor Immune Estimation Resource (TIMER) to analyze and visualize the association between immune infiltrates and expression level of the identified potent antigens. Consistency matrix analysis was used to identify corresponding gene modules and immune subtypes. Univariate and multivariate Cox regression as well as log-rank tests were used to assess the prognostic distinctions of the GBM immune subtypes. We used the reduce Dimension function of the Monocle package to visualize the distribution of immune subtypes across TCGA and CGGA patients.

Results

Different cancer have different potential RNA Antigens. In GBM, seven tumor antigens (DIRAS3, EGFLAM, FUCA1, LILRB2, MPZL2, OSMR, and STC1) were identified as promising candidates for developing anti-GBM mRNA vaccines. Based on the immune gene expression profiles, GBM was classified into four immune subtypes (IS1–IS4). IS1 and IS4 were immunologically “hot,” which means they are more likely to benefit from vaccines. In addition, six immune hub genes (CLCF1, CD54, SOCS3, PLAUR, LIF, and BCL3) were identified as potential biomarkers for GBM vaccine effectiveness.

Conclusion

Different cancer types should have different potential RNA vaccine. For instance, DIRAS3, EGFLAM, FUCA1, LILRB2, MPZL2, OSMR and STC1 were identified to be related with poor survival and antigen presenting cell (APC) infiltration in GBM, and they are promising targets for mRNA vaccine development. GBM patients could be divided in to four immune subtypes based on the clustering of immune-related gene expression, and those in the IS1 and IS4 groups were candidates for vaccination.

GEPIA analysis

As described by Huang et al.

^14,

¹⁵

, we used Gene Expression Profiling Interactive Analysis (GEPIA) version 2 (GEPIA2, https://gepia2.cancer-pku.cn), which is an open-access online tool that collects tumor and normal samples from The Cancer Genome Atlas (TCGA, (https://www.cancer.gov/tcga) and the Genotype-Tissue Expression databases, to perform the differential expression analysis of GBM genes and to calculate the prognostic index of each selected antigen. The GEPIA2 web server also provides extended gene expression quantification from the gene level to the transcript level and supports analysis of a specific cancer subtype as well as comparison between subtypes

¹⁹

. In this study, we used analysis of variance to identify the significantly differentially expressed genes based on |log2FC| values > 1 and q values < 0.01. We separated the patients with a 50% (median) expression level cutoff and performed overall survival (OS) and relapse-free survival (RFS) analysis using the Kaplan-Meier method. We used the Cox proportional hazards regression model to calculate the hazard ratio. The parameter setting was consistent in each analysis without adjustment for any P value. We considered P values < 0.05 to be statistically significant.

cBioPortal analysis

We used the cBio Cancer Genomics Portal (cBioPortal, https://cbioportal.org) to visualize the gene alteration of potential antigens against GBM in the TCGA. cBioPortal allows interactive exploration of multidimensional cancer genomics data sets

²⁰

. Statistical significance was set at P < 0.05.

TIMER analysis

The Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) is an open-access resource for comprehensively investigate molecular characterization of tumor-immune interactions

²¹

. In this study, we used TIMER to analyze and visualize the association between immune infiltrates and expression level of the potent antigens we identified. P < 0.05 was considered to be statistically significant.

Data extraction of immune-related gene

The normalized gene expression data from 325 patients and clinical data from 306 patients were collected from the Chinese Glioma Genome Atlas (CGGA, https://www.cgga.org.cn/tools.jsp), and RNA-seq data from 167 patients and clinical data from 600 patients were collected from the TCGA. We extracted a total of 2006 immune-related genes from both the discovery and validation cohorts, which included single immune cell-specific, co-stimulating, and co-suppressing molecule-related genes, cytokines and cytokine receptor-related genes, antigen processing and presenting-related genes, and other immune-related genes. TCGA and CGGA tumor samples lacking clinical information and normal tissue data were excluded, as were genes with 0 transcripts per million (TPM) in > 50% of the samples. Ultimately, 103 immune cell-related genes with log2 (TPM+1) fold-change were included in the subsequent analysis.

Discovery and validation of the immune subtypes

Immune-related genes were clustered on the basis of their expression profiles, and we developed a consistency matrix to identify corresponding gene modules and immune subtypes. We applied the partition around medoids algorithm using the “1-Pearson correlation” distance metric and performed 500 bootstraps, each involving 80% of the patients in the discovery cohort. We set cluster varied from 2 to 8 to determine the optimal number of clusters, and we defined the optimal partition by evaluating the consensus matrix and the consensus cumulative distribution function. We used an independent TCGA cohort to validate the immune subtypes with the same settings. The consistency of immune subtypes between discovery and validation datasets was quantified by calculating the Pearson correlation and the proportion of the intra-group in the centroid of the gene module scores.

Prognostic evaluation of immune subtypes and construction of immune landscape

The prognostic distinctions of the GBM immune subtypes were assessed using univariate and multivariate Cox regression and log-rank tests with stage and grade as covariates. We used analysis of variance to assess the relationship between immune subtypes with different immune-related molecular and cellular characteristics. Gene Ontology analysis was performed using The Database for Annotation, Visualization and Integrated Discovery v6.8 (https://david.ncifcrf.gov/) with the genes expressed (TPM >1) as the background

²²

. We used single-sample gene set enrichment analysis (ssGSEA), which is the measure of genes that are coordinately up- or down-regulated within a sample, to calculate immune enrichment scores for each sample. We screened the most frequently mutated genes using the chi-square test.

We performed dimensionality reduction analysis to visualize the distribution of immune subtypes across TCGA and CGGA patients using the reduce Dimension function of the Monocle package with a Gaussian distribution. According to the methods published earlier by Huang et al.

^14,

¹⁵

, the maximum number of components was set to 4. Finally, we used the function plot cell trajectory with color-coded immune subtypes in the Monocle package to visualize the immune landscape. Weighted correlation network analysis was used to constructed co-expression modules of immune genes according to the methods published earlier by Huang et al.

^14,

¹⁵

.

Identification of potential tumor antigens of different cancer type

Using TCGA data, we analyzed the expression levels of natural receptors and ligands across different tumors and found that different tumors possess distinct natural receptors and ligands (Fig. 1). This indicates that different tumors should have different RNA vaccines. Among these, we found that gliomas have unique receptor and ligand expression. Therefore, we will focus our subsequent research on the natural receptors and ligands of gliomas.

Identification of potential tumor antigens of GBM

We screened the aberrantly expressed genes of GBM to identify potential antigens. Among the 7659 differentially expressed genes in GBM, we detected 5221 overexpressed genes that may encode tumor-associated antigens (Fig. 2A). Subsequently, 8182 mutated genes encoding tumor-specific antigens were screened by analyzing altered genome fraction and mutation counts in individual patients (Fig. 2B, C). Genes with the highest alteration frequency in the fraction genome-altered group were ASAH1, DES, DIS3, FAM86B2, PROSER1, TUBA4A, EGFR, and SEC61G (Fig. 2D). High mutation counts were observed in CDKN2A-DT, SEC61G, MTAP, CDKN2A, CDKN2B, LANCL2, EGFR, VOPP1, PTEN, and DMD (Fig. 2E). Overall, we identified 1633 overexpressed and frequently mutated tumor-specific genes.

Identification of tumor antigens associated with GBM prognosis and antigen-presenting cells

We next screened the aforementioned genes for prognosis-related tumor antigens that may have immune stimulatory or inhibitory effects as candidates for developing mRNA vaccines. Among the 1633 overexpressed and frequently mutated tumor-specific genes, 35 genes were closely associated with the OS of GBM patients, and 7 of those genes showed significant correlation with RFS (Fig. 3A). Patients overexpressing the seven genes (DIRAS family GTPase 3 (DIRAS3), alpha-L-fucosidase 1 (FUCA1), leukocyte immunoglobulin like receptor B2 (LILRB2), myelin protein zero like 2 (MPZL2), oncostatin M receptor (OSMR), stanniocalcin 1 (STC1), and EGF like, fibronectin type III and laminin G domains (EGFLAM)) in the tumor tissues had significantly shorter survival time compared to the group with lower expression (Fig. 3B–H). Thus, these seven gene candidates were identified as being critical for GBM development and progression. What’s more, higher expression levels of DIRAS3, FUCA1, and LILRB2 were significantly associated with increased tumor infiltration of B cells, macrophages, and dendritic cells (Fig. 4A, C, D). EGFLAM, MPZL2, OSMR, and STC1 expression levels also tended to be associated with increased infiltration of immune cells, but the patterns were more variable (Fig. 4B, 3E–G). All in all, these findings suggest that the seven identified tumor antigens may be directly processed and presented by antigen presenting cells (APCs) to T cells and recognized by B cells to trigger an immune response. Therefore, they are promising candidates for developing mRNA vaccines against GBM.

Identification of potential immune subtypes of GBM

Immunophenotyping is helpful for identifying appropriate patients for vaccination since it can reflect the immune status of the tumor and its microenvironment. So, we analyzed the expression profiles of 402 immune-related genes in 325 GBM samples from the CGGA database to construct consensus clustering. We chose k = 4, where immune-related genes appeared to be stably clustered according to their cumulative distribution function and their function delta area (Fig. 5A, B), and we obtained four immune subtypes designated as IS1–IS4 (Fig. 5C). IS3 was associated with the best prognosis, whereas IS1 and IS2 had poor survival probability (Fig. 5D). We also analyzed the expression profiles of 396 immune-related genes in 167 GBM samples from the TCGA database to construct consensus clustering. Consistent with the results obtained with the CGGA cohort, four immune subtypes were obtained and related to prognosis in this cohort (Fig. 5E).

We investigated the subtype distribution in samples with different 1p19q codeletion status and MGMTp methylation status in the CGGA cohort and found that IS2 was only distributed in those with 1p19q non-code status (Fig. 5F), whereas IS1–IS4 were irregularly clustered in samples with both MGMTp methylated and un-methylated status (Fig. 5G). In the TCGA cohort, pan-CDKN2A and TP53 were irregularly distributed in IS1–IS4 (Fig. 5H). Furthermore, we investigated the IS1–IS4 distribution across classical, mesenchymal, neural, and proneural subtypes, which are used for gene expression-based molecular classification of GBM. The classical subtype was only associated with IS2, the mesenchymal subtype was only associated with IS1, and the proneural subtype was only associated with IS3 (Fig. 5I). These results show that immunotyping can be used to predict prognosis of GBM patients and classical subtypes.

Association of immune subtypes with tumor mutational burden and mutational status

Higher tumor mutation burden (TMB) is correlated with stronger anti-cancer immunity. We calculated the TMB in each patient using the mutect2-processed mutation dataset of TCGA. However, there was no significant difference in TMB among the four immune subtypes (Fig. 6A). Eight genes (TP53, EGFR, PTEN, TNN, SPTA1, RYR2, MUC16, and IDH1) were most frequently mutated in these subtypes (Fig. 6B). These findings suggest that the tumor antigens do not differ significantly among the IS1–IS4 subtypes.

Association between immune subtypes of GBM and immune modulators

Both immune checkpoints (ICPs) and immunogenic cell death (ICD) modulators play critical roles in modulating host anti-tumor immunity, which could influence the efficacy of mRNA vaccines. Therefore, we analyzed their expression levels in the different subtypes. Forty-three ICP-related genes were detected in both the CGGA and TCGA cohorts, of which 40 (93%) genes in the CGGA cohort (Fig. 7A) and 39 (91%) genes in the TCGA cohort (Fig. 7B) were differentially expressed among the four immune subtypes. Furthermore, the overall expression level of ICPs in the TCGA cohort was higher than that in the CGGA cohort. Likewise, 21 ICD genes were detected in both the CGGA and TCGA cohorts, of which 20 (95%) genes in the CGGA cohort (Fig. 7C) and 19 (90%) genes in the TCGA cohort (Fig. 7D) showed significant differences among the four immune subtypes. Therefore, immunotyping can reflect the expression levels of ICPs and ICD modulators, which may prove to be useful therapeutic biomarkers of the effectiveness of vaccines.

Association between immune subtypes and GBM-related tumor marker

Neuron-specific enolase (NSE) is the most commonly used prognostic tumor biomarker for neuroma, and a high value indicates cancer progression, poor prognosis, or cancer relapse

²³

. In this study, both the CGGA and TCGA cohorts displayed significant differences in NSE expression levels across the immune subtypes (Fig. 8A, B). For example, IS4 showed lower NSE expression in the CGGA cohort (Fig. 8A), whereas IS3 showed lower NSE expression in the TCGA cohort (Fig. 8B). These outcomes demonstrate that immune subtype is a better indicator than NSE for anticipating GBM patient prognosis.

Cellular and molecular characteristics of immune subtypes

Because the tumor immune status to a great extent will impact the efficacy of mRNA vaccines, we further described the resistant cell segments in the four immune subtypes by scoring 28 marker genes in both the CGGA and TCGA cohorts utilizing ssGSEA. The immune cell composition differed significantly among the subtypes in the CGGA cohort (Fig. 9A). Most immune cells were more enriched in the IS1 and IS4 immune subtypes. For instance, the scores of activated CD8 T cells and immature B cells were significantly higher in IS1 and IS4 compared to IS2 and IS3 (Fig. 9A, C). Therefore, IS1 and IS4 are immunologically “hot” and IS2 and IS3 are immunologically “cold” phenotypes. Similar trends were seen in the TCGA cohort, as most immune cells were more enriched in the IS1 and IS4 immune subtypes (Fig. 9B, D). These results suggest that immune subtype reflects the GBM immune status and that it can be used to identify suitable patients for mRNA vaccination. The mRNA vaccines with these antigens should induce immune infiltration in patients with immunologically “cold” (IS2 and IS3) tumors.

Immune landscape of GBM

The immune gene expression profiles of each patient were used to develop the immune landscape of GBM for mRNA vaccine application (Fig. 10A). The horizontal axis was correlated with various immune cells, of which effector memory CD8 T cells, type 1 T helper cells, natural killer T cells, mast cells, gamma delta T cells, central memory CD4 T cells, and activated CD8 T cell showed the most negative correlation, whereas the vertical coordinate was mostly negatively associated with central memory CD8 T cells (Fig. 10B). IS1, IS2, and IS4 were further divided into two subgroups according to the distribution location of the three immune subtypes in the immune landscape (Fig. 10C), and several immune cells differed significantly among these subsets (Fig. 10D). For instance, the enrichment scores of activated B cells, activated CD8 T cells, immature B cells, macrophages, MDSCs, natural killer cells, and T follicular helper cells in IS1B were lower than those in IS1A. Likewise, IS2A scored lower than IS2B in terms of activated B cells, activated CD4 T cells, activated CD8 T cells, CD56bright natural killer cells, CD56dim natural killer cells, effector memory CD4 T cells, memory B cells, monocytes, plasmacytoid dendritic cells, and type 2 T helper cells. IS4B scored lower in terms of effector memory CD8 T cells, macrophages, and MDSCs, while IS4A scored lower in terms of monocytes and type 2 T helper cells. Thus, mRNA vaccines may be relatively viable and more effective in IS1B and IS2A. Furthermore, patients were separated into six groups based on their location in the immune landscape and were subjected to prognostic comparison (Fig. 10E). Patients in group 5 showed the poorest survival probability (Fig. 10F), which is consistent with the aforementioned results. Overall, these findings suggest that the immune landscape based on immune subtypes can be used to define the immune components of each GBM patient and predict prognoses that help in selecting suitable patients and personalized therapeutics for mRNA vaccines.

Identification of immune gene co-expression modules and hub genes of GBM

Co-expression modules of immune genes were constructed by clustering the samples using weighted correlation network analysis (Fig. 11A). We set the soft threshold at 4 for the scale-free network (Fig. 11B). We calculated eigengenes of each module and merged the close modules into a new one (Fig. 11C). We further analyzed the distribution of the eigengenes in the four immune subtypes of these modules and detected significantly different distributions in every module (Fig. 11D). IS3 had the lowest number of eigengenes in the black, blue, pink, brown, and red modules, and IS2 had the lowest number of eigengenes in the magenta, turquoise, yellow, and green modules. IS1 had the highest number of eigengenes in the black, blue, pink, and turquoise modules, and IS4 had the highest number of eigengenes in the magenta, brown, red, and yellow modules. Thus, IS3 and IS2 corresponded to immunologically cold and IS1 and IS4 to inflamed tumors.

Further prognostic correlation analysis showed that the blue module was significantly associated with the prognosis of GBM (Fig. 12A). Moreover, genes in the blue module were enriched with cytokine-cytokine receptor interaction, chemotaxis, interleukin-related, and immune-related pathways and biological processes (Fig. 12B). The blue module also was negatively related to component 1 of the immune landscape (Fig. 12C). Analysis of the prognostically relevant genes of the blue module showed that higher expression scores correlated with better prognosis in the CGGA and TCGA cohorts, which is consistent with the abovementioned prognosis findings (Fig. 12D). The immune-related genes CLCF1, CD54, SOCS3, PLAUR, LIF, and BCL3 were associated with the blue module, and we conclude that these six immune hub genes are potential biomarkers for effectiveness of mRNA vaccines in GBM.