Prognostic prediction and immune landscape analysis based on m6A methylation modification and vasculogenic mimicry in cervical cancer

2.1 Data acquisition

RNA sequencing data (tragments per kilobase of transcript per million fragments mapped (FPKM) values) were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database, which included 306 tumor samples and three normal samples. The FPKM values were first converted to transcripts per kilobase million (TPM) values. Next, we downloaded the corresponding clinical information while excluding samples with a survival time of 0. Finally, 294 patients were included in our analysis. We also downloaded copy number variation data and acquired somatic mutation data from UCSC Xena (https://xenabrowser.net).

2.2 Detection of mVM genes

Twenty-eight MAGs and 25 VRGs were collated from previous literature, thus yielding a total of 53 mVM genes [24, 25, 26, 27]. Pearson correlation analysis was then used to assess the correlation between MAGs and VRGs. Gene interactions were analyzed by the online STRING tool (https://cn.string-db.org/) and Cytoscape software version 3.9.0 was used to visualize the correlation results.

2.3 Unsupervised clustering analysis

When investigating the specific role of mVM genes in CC, we determined the optimal number of clusters to cluster patients into different subgroups using the “ConsensusClusterPlus” tool in the R package [28].

2.4 Construction of the mVM prognostic signature

CC patients were randomized into training and validation sets by the “caret” package and used univariate Cox analysis to screen the prognosis-related genes. Subsequent LASSO regression led to the identification of eight key genes thus allowing us to establish a mVM-related prognosis signature. Patients were divided into high- and low-risk groups based on the median risk score (referred to as the mVMscore), and the difference in prognosis between the two groups was assessed using Kaplan-Meier analysis. In addition, time-dependent curves were plotted to evaluate the accuracy of the signature at different periods. The mVMscore, age, grade and TNM (tumor, node, metastasis) staging metrics were incorporated into ROC curves to further evaluare the prognostic predictive ability between mVMscores and clinical factors.

2.5 Functional enrichment analysis and gene set variation analysis (GSVA)

The “clusterProfiler” package was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to identify functional and pathway enrichment [29]. The “c2.cp.kegg.v7.5.symbols” gene set was downloaded from the Molecular Signatures Database (MSigDB), and differences in biological processes under different mVM modification patterns were determined by GSVA [30].

2.6 Tumor microenvironment and immune infiltration analysis

To quantify tumor-infiltrating immune cells, we used the CIBERSORT algorithm [31] to determine the proportion of infiltrating immune cell types. Single sample gene set enrichment analysis (ssGSEA) was then used to evaluate the abundance of infiltrating immune cells. The “ESTIMATE” tool in the R package was then used to calculate immune scores, stomal scores, estimate scores and tumor purity in patients with CC [32]. Tumor immune dysfunction and rejection (TIDE) scores were then downloaded from the TIDE database (http://tide.dfci.harvard.edu/) and used to evaluate the likelihood of a response to immunotherapy.

2.7 Statistical analysis

All statistical tests were performed in R software version 4.3.1. Pearson’s correlation analysis was used to investigate the association between MAGs and VRGs. Kaplan-Meier survival analysis, with the log-rank test, was then used to compare patients in different subgroups. Statistical significance was set at p < 0.05.

3.1 Multi-omics characterization of the landscape of m6A regulators in CC

To reveal the regulatory role of m6A modifications in CC, we first explored the genetic changes of 28 m6A regulators. Our analysis showed that the overall mutation frequency of m6A regulators was low, with somatic mutations occurring in 50 of the 289 samples (17.3%); the highest frequencies of mutation (3%) were detected in Leucine-rich pentatricopeptide repeat containing protein (LRPPRC), ZC3H13 and YTHDC2 (Fig. 1A). Fig. 1B depicts the majority of m6A regulator copy number deletions, with only Fragile X mental-retardation protein (FMR1), RBMX, HNRNPC, METTL3, YTHDC1 and VIRMA exhibiting amplification. Fig. 1C shows the location of copy number variations for m6A regulators on different chromosomes. In addition, to evaluate the interaction between genetic variations and the expression of m6A regulators, we further analyzed differences in the expression levels of m6A regulators between tumor and normal samples. We found that METTL16, FTO, RBM15, YTHDF2 and HNRNPA2B1 exhibited differential expression when compared with normal tissues; specifically, METTL16 and FTO exhibited reduced expression levels in tumor tissues, while RBM15, YTHDF2 and HNRNPA2B1 exhibited elevated expression levels (Fig. 1D). Pearson correlation analysis confirmed the close correlation between m6A regulators (Fig. 1E). Fig. 1F demonstrates the prognostic value of the three RNA modification species in patients with CC, thus forming a complex regulatory network. Collectively, these results revealed that m6A regulators may play an important role in the progression of CC.

Fig. 1.Multi-omics characterization landscape of m6A regulators in cervical cancer. (A) The mutation landscape of 28 m6A regulators. (B) Copy number variation alteration frequency of m6A regulators. (C) The location of CNV alteration of m6A regulators on chromosomes. (D) Differential expression of m6A regulators in tumor and normal samples. *p < 0.05, **p < 0.01. (E) Pearson correlation analysis reveals positive and negative correlations among m6A regulators. (F) Prognostic network diagram of the three RNA modification species. The left half of circle in different colors represents different RNA modification types. Lines between m6A regulators represent their positive/negative correlations. The size of the circles represents the relationship between m6A regulators and survival. The larger the circle the more likely it is to be a prognosis-related gene.

3.2 Association between m6A modifications and VM

A large body of evidence suggests that m6A modifications are involved in regulating the progression of a variety of cancers, and that VM is closely associated with tumor invasion and migration. To investigate whether there is a correlation between m6A modifications and VM, we confirmed the association within VRGs by Pearson’s correlation analysis; this analysis indicated that Rho kinase type 1 (ROCK1) exhibited the strongest correlation with Mitogen-activated protein kinase 1 (MAPK1), with a correlation coefficient of 0.84 (Fig. 2A). In addition, we generated a heatmap which illustrated significant correlations between m6A modifiers and VRGs (Fig. 2B). We also generated a protein-protein interaction (PPI) network which further confirmed the intricate interactions between m6A regulators and VRGs (Fig. 2C).

Fig. 2.Correlation between m6A regulators and VM. (A) The expression correlation analysis of VRGs in CC. (B) The correlations between m6A regulators and VRGs. *p < 0.05, **p < 0.01, ***p < 0.001. (C) PPI network revealed complex relationships between MARs and VRGs.

3.3 Identification of mVM molecular patterns and biological function analysis

First, we performed univariate Cox regression analysis on the 53 mVM genes (Supplementary Fig. 1) and selected significant genes (p < 0.05) for unsupervised consensus clustering. We found that k = 3 was the most appropriate choice when cluster stability increased from k = 2 to 9 (Fig. 3A–C); this led to the classification of patients into three subgroups: cluster C1, cluster C2 and cluster C3. When considering the three subgroups, survival analysis showed that cluster C1 had a significantly shorter survival and the worst prognosis (Fig. 3D).

Next, we used GSVA to further compare biological behaviors between the three subgroups. cluster C1 was predominantly enriched in focal adhesion, glycosaminoglycan degradation, glycosaminoglycan biosynthesis, chondroitin sulfate, galactose metabolism and the nod-like receptor signaling pathway. cluster C2 was significantly enriched in the linoleic acid metabolism pathway while cluster C3 was associated with metabolic or biosynthetic activation phenotypes (Supplementary Fig. 2A–C). Analysis of immune cell infiltration further showed that cluster C1 was mainly enriched in innate immune cells, such as γδ T cells, mast cells, neutrophils, regulatory T cells and Type 2 T helper cells. Activated B cells accounted for the largest proportion of cells in cluster C2, while cluster C3 featured the highest number of immature dendritic cells (Fig. 3E).

Fig. 3.Recognition of mVM genes modification patterns. (A) Consensus clustering matrix at k = 3. (B) Cumulative Distribution Function (CDF) curve for consensus clustering analysis at k = 2–9. (C) Relative change in area under the CDF curve at k = 2–9. (D) Survival analysis of different clusters by Kaplan-Meier. (E) The abundance of immune-infiltrating cells in the three clusters. *p < 0.05, **p < 0.01, ***p < 0.001.

3.4 Construction and validation of prognostic signature

With the purpose of better predicting the clinical prognosis based on mVM genes, patients with CC in the TCGA database were randomly divided into a training set and a validation set for model construction and validation. First, we performed univariate Cox regression analysis on 53 mVM genes (Supplementary Fig. 1) and identified 14 prognosis-related genes (p < 0.05). To improve the predictive effect, we used the LASSO-Cox regression algorithm for further screening; this strategy ultimately identified eight prognostic genes (Fig. 4A,B). Kaplan-Meier analysis demonstrated that high expression levels of all eight genes identified by LASSO analysis were significantly associated with a poor patient prognosis (Supplementary Fig. 3). In addition, we constructed a prognostic signature by calculating the risk score for each patient, which we defined as the mVMscore. Patients with CC were then categorized into low-risk and high-risk groups based on the median mVMscore. Fig. 4C–E shows the mVMscore, survival time, live status and a heatmap of mVM genes expression in CC. Kaplan-Meier survival curves further showed that the high mVMscore group had a worse prognosis than the low mVMscore group, both in the training group (Supplementary Fig. 4A), the validation group (Supplementary Fig. 4B), and in the entire cohort (Fig. 5A).

Fig. 4.Construction of a prognosis signature base on mVM genes. (A,B) LASSO Cox analysis identified eight prognosis-related genes. (C–E) Distribution of mVmscore, survival status of CC patients and heatmap expression of characteristic gene expression in the training cohort (C), testing cohort (D) and entire TCGA cohort (E).

Then, in order to better understand whether the mVMscore could effectively predict the prognosis of a given patient, we plotted time-dependent ROC curves based on mVMscore. This allowed us to determine the change of survival model and prediction accuracy at different time points by calculating the areas under the curve (AUCs) at 1 year, 3 years and 5 years. Analysis showed that the 1-year, 3-year and 5-year area under the curve (AUC) values in the training set were 0.713, 0.735 and 0.706, respectively (Supplementary Fig. 4C). The AUC values for 1-, 3- and 5-years were all above 0.7 for the testing cohort (Supplementary Fig. 4D) and entire TCGA dataset (Fig. 5B). In addition, the 1-, 3- (Supplementary Fig. 4E,F) and 5-year ROC curves (Fig. 5C) for all cervical cancer patients showed that the prognostic ability of the mVMscore was significantly better than other clinical traits, such as grade, age and TNM stage. Therefore, our new prognostic model has good prognostic value.

3.5 Analysis of the correlation between mVMscore, clinical traits and somatic mutations

The Sankey diagrams demonstrated the analysis workflow of mVMscore, cluster and survival status (Fig. 5D). We observed that all three subgroups were correlated with patient survival status, mVMscore was higher in cluster C1; both cluster C2 and cluster C3 were correlated with a low mVMscore (Fig. 5E). We also observed that mVMscore was higher in patients with death status than in patients with survival status (Fig. 5F). Subsequently, we investigated the association between Tumour mutation burden (TMB) and mVMscore and found that patients in the high TMB group had a significantly better prognosis (Fig. 5G). When the high TMB group was combined with the low mVMscore group, we found that survival was significantly better than for the other groups (Fig. 5H).

Fig. 5.Analysis of mVMscore correlation with clinical traits and somatic mutations. (A) Kaplan-Meier curves of the mVMscore in the entire TCGA dataset. (B) Time-dependent ROC curve of the mVMscore prognostic signature for 1-, 3- and 5-year overall survival (OS) of CC patients in the entire TCGA dataset. (C) ROC curve analysis of the mVMscore and clinicopathological parameters (grade, age, T stage, N stage and M stage) for 5-year OS. (D) The alluvial diagram depicted the changes of clusters, mVMscore and survival status. (E) Comparison of mVMscore among different clusters. (F) Comparison of mVMscore between different survival states. (G) Kaplan-Meier curves for high-TMB vs. low-TMB groups. (H) Kaplan-Meier survival analysis for patients with different expression levels of mVMscore & TMB.

Next, we compared somatic mutation frequencies between the high mVMscore and low mVMscore subgroups and found that overall mutation frequencies did not differ significantly between these two groups. The most common type of mutation was nonsense mutations; furthermore, the somatic mutation rates of TTN (36% vs. 22%) and Phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha (PIK3CA) (30% vs. 23%) were significantly higher in the high mVMscore subgroup (Fig. 6A) than in the low mVMscore subgroup (Fig. 6B). In addition, we also demonstrated that the top 20 mutated genes were significantly less frequent in the high mVMscore group (Fig. 6C), but were more frequent in the low mVMscore group (Fig. 6D).

Fig. 6.Comparison of somatic mutations among different subgroups. (A) The mutation alteration in high-risk group. (B) The mutation alteration in low-risk group. (C) Correlation of mutated genes in high mVMscore groups. (D) Correlation of mutated genes in low mVMscore groups.

3.6 Correlation between mVMscore and the tumor immune landscape

To enhance our understanding of the tumor microenvironment, we next assessed the abundance of immune-infiltrating cells in the tumor microenvironment by performing CIBERSORT and ssGSEA analyses. These analysis showed that immune cells were significantly enriched in the high-risk group, including activated B cells, activated cytotoxic CD8 T cells, eosinophils, immature B cells, immature dendritic cells, myeloid-derived suppressor cells (MDSCs), macrophages, regulated immature B cells, regulatory T (Treg) cells, follicular helper T (Tfh) cells, type 1 helper T (Th1) cells and natural killer (NK) cells (Fig. 7A,B). Based on the ESTIMATE algorithm, we then analyzed the differences between stromal and immune cells in the high and low-mVMscore groups. Analysis showed that the immune score, estimated score, and stromal score were higher in the low-risk group than in the high-risk group (Fig. 7C). TIDE scores were higher in the low-mVMscore group (Fig. 7D), thus indicating that immune escape was more likely to occur in patients in the low-risk group, and that the efficacy of immunotherapy was greatly reduced. Therefore, we also analyzed the expression levels of immune checkpoint genes, including Programmed death-ligand 1 (PD-L1), Programmed cell death protein 1 (PD-1) and CTL-associated antigen 4 (CTLA4). Analysis showed that the expression levels of immune checkpoint genes were up-regulated in all of the low-mVMscore groups (Fig. 7E). These findings suggested that immunotherapy efficacy was weaker in the low-risk population, and that the benefit of immunotherapy may be greater in the high-mVMscore group. Collectively, our results suggest that the mVMscore has great potential for assessing the tumor microenvironment, immune function, and immunotherapy in patients with CC.

Fig. 7.Potential for predicting immunotherapeutic response in mVM signature. (A) Differences in infiltration of immune cell types between the two subgroups. (B) The infiltrating levels of different immune cell types. (C) ESTIMATE was performed to compare immune scores, stromal scores, tumor purity and estimate scores in two groups. (D) TIDE prediction difference in the high- and low-risk patients. (E) Violin plot depicted the expression of checkpoint genes in the high- and low-mVMscore groups. *p < 0.05, **p < 0.01, ***p < 0.001.