The transcriptome profiling data of patients with cervical cancer that we analyzed were retrieved from the TCGA database (https://portal.gdc.cancer.gov/) on May 30, 2022. The data included 309 RNA-seq files. The tumor-corresponding clinical characteristics of cervical cancer were also downloaded from TCGA, and the data included 304 cases of tumors. The GSE52903 cohort was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which includes 55 cervical cancer RNA-seq files and corresponding clinical characteristics.

All of the transcriptome data were analyzed by using the “ESTIMATE” R package (https://www.r-project.org, v1.0.13) via Rstudio (https://www.rstudio.com, v2022.02.2+485). The calculated immune, stromal, and ESTIMATE scores were observed to display a component ratio in TME immune cell infiltration. ESTIMATES scores consisted of both ImmuneScores and StromalScores. Higher scores indicated a larger number of components in the TME.

The samples were divided into the “high scores” group and “low scores” group according to the medians of the ImmuneScores, StromalScores, and ESTIMATEScores. A Kaplan–Meier survival analysis was performed in two groups via the R package “survival” (https://www.bioconductor.org/, v3.3-1) and “survminer” (https://www.bioconductor.org/, v0.4.9). The primary prognostic endpoint was overall survival (OS). In addition, some clinical characteristic data were permitted to assess the associations with the scores, which included age, stage (defined by FIGO), histopathological grades, and Tumour, Nodes, and Metastasis (TNM) stage. These comparisons were carried out using the Kruskal-Wallis rank-sum test, with p 0.05 considered statistically significant.

DEGs were distinguished by using the “high scores” group and “low scores” group via the R package “limma” (https://www.bioconductor.org/, v3.48.3). A gene was only regarded as a DEG if two conditions were met: (1) false discovery rate (FDR) 0.05 and (2) $|$log2(fold change)$|$ $>$1. The R packages “VennDiagram” (https://www.bioconductor.org/, v1.7.3) and “pheatmap” (https://www.bioconductor.org/, v1.0.12) were used for DEG visualization.

Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis with 791 DEGs were performed by using the R package “clusterProfiler” (https://www.bioconductor.org/, v4.0.5) and visualized by using the website HIPLOT (https://hiplot.com.cn) and GOBubble Plot module (https://www.bioconductor.org/, v0.1.0).

The protein and protein $>$0.4 and homo protein interaction (PPI) networks were constructed on the website STRING (https://cn.string-db.org/). All of the DEGs were screened with confidence scores, and interactions based only on experimental evidence were selected as the organism to build the network. The network was visualized using Cytoscape (https://cytoscape.org/, v3.8.0). Moreover, the modular analysis was performed using the Cytoscape tool cytoHubba (https://www.bioconductor.org/, v2.5.3), and the top 10 genes with at least 15 adjacent nodes were selected according to the multinetwork clustering method (PPI key genes). We also performed a COX regression analysis with the R package “survival” to explore the relationship between each DEG and OS. Similarly, p 0.05 was considered the cutoff criteria, wherein we obtained 147 genes related to patient survival status (the COX key genes). Subsequently, the R package “VennDiagram” was utilized to portray the intersection of key COX and PPI genes.

The concerned gene GRAP2 was selected. Afterward, we summarized the correlation between clinical information and GRAP2 expression of patients in the TCGA database, Pearson’s Chi-squared test, Fisher’s exact test, and Wilcoxon rank sum test were used to test the p-value. The expression of each gene was matched to clinical features using “ggpubr” (https://www.bioconductor.org/, v0.4.0) and “survival”. For these comparisons, a Wilcoxon rank-sum was employed, and p 0.05 was considered statistically significant. To validate the reliability of the results, we extracted the GSE52903 data of cervical cancer patients from the GEO database. Due to the fewer number of samples (only 55 samples), we used the R package “maxstat” (https://www.bioconductor.org/, v0.7-25) to calculate the best cutoff value of $|$log2(fold change)$|$, wherein the minimum grouping sample size was set greater than 25%, and the maximum sample size grouping was less than 75%. Based on the best cutoff value of 5.81137, the patients were divided into two high ($>$5.81137) and low ($\le$5.81137) groups. By using the R package “survival” we performed a Kaplan–Meier survival analysis to estimate survival differences between the two groups.

To further explore the role of GRAP2 in cervical cancer, GSEA, a computational method that determines whether an a priori defined set of genes shows statistically, was performed by using GSEA ( http://www.gsea-msigdb.org/gsea/index.jsp, v4.2.3) software to analyze the gene matrix of GO and KEGG associated with the gene GRAP2 expression level. Moreover, p 0.05 was considered to be significantly enriched.

The website “CIBERSORx” (https://cibersortx.stanford.edu/) provides the CIBERSOR tool to evaluate the content of immune cells in tumor cases. The Wilcoxon rank-sum test and Pearson’s correlation were used to assess the proportions and gene expression. The results were visualized by using the “vioplot” (https://www.bioconductor.org/, v0.3.7), “ggplot2” (https://www.bioconductor.org/, v3.3.2), and “ggpubr” R packages.

The current study was performed by using the following analysis process (Fig. 1). First, we downloaded the RNA-seq profiles of 304 cervical cancer cases and the corresponding clinical information from the TCGA and GEO databases. Afterward, calculations of TIC ratios and immune and stromal proportions were carried out separately using the CIBERSORT and ESTIMATE algorithms. We identified DEGs between groups with low and high ImmuneScore and StromalScore. Simultaneously, clinical information and TME scores were integrated for univariate Cox regression analyses, and we found that 147 genes showed a significant association with prognosis (the Cox key genes). Subsequently, the DEGs were subjected to PPI network analysis and univariate Cox regression analyses, and 30 genes were selected by the cytoHubba plugin (the PPI key genes). Finally, an intersection analysis of Cox key genes and PPI key genes was used to obtain 12 prognosis key genes. The gene GRAP2 attracted our attention; therefore, GRAP2 was further analyzed, including expression, survival, and clinicopathological characteristics, as well as GSEA and correlations with TICs proportions.

Fig. 1.

Overview of this study’s analytical flowchart.

Cervical cancer patients’ immune infiltration scores were evaluated for prognoses (Fig. 2). For each of the TME scores (ImmuneScore, StromalScore, and ESTIMATEScore), half of the 304 tumor samples (with complete clinical data) were classified into the high score group (152 cases), and the other half of the samples were classified into the low score group (152 cases) based on each median score. Thereafter, Kaplan–Meier survival curves were constructed.

Fig. 2.

A comparison of ImmuneScore, StromalScore, and ESTIMATEScore with survival among TCGA-CESC participants. (A) ImmuneScore (p = 0.072, HR = 0.65, 95% CI = 0.40–1.04), (B) StromalScore (p = 0.484, HR = 0.85, 95% CI = 0.53–1.35), and (C) ESTIMATEScore (p = 0.024, HR = 0.0.58, 95% CI = 0.36–0.94) survival analyses for cervical cancer patients with low and high scores. HR, hazard ratio; CI, confidence interval.

OS was better for patients with a high ESTIMATEScore (log-rank test p = 0.024, hazard ratio (HR) = 0.58) than the low scores (Fig. 2C). However, the StromalScore and ImmuneScore alone were not significantly associated with OS (Fig. 2A,B). According to these findings, the entire cervical TME might play an important role in cervical cancer prognosis.

We also analyzed the relationship between TME scores and clinical characteristics (Fig. 3). The latter parameter included age, FIGO stage, pathologic grade, and the TNM stage. Along with the progression of the M classification, the StromalScore and ESTIMATEScore notably declined (except for the ImmuneScore) (Fig. 3N,O). This result is different from that of a previous report [16] due to more samples being included and more cases being assessed by our team. These findings clarified that the TME of cervical cancer plays a crucial role in the invasion and metastasis of cancer.

Fig. 3.

Associations of ImmuneScore, StromalScore, and ESTIMATEScore with age, FIGO stage, grade, and TNM staging. (A–C) Kruskal-Wallis rank sum test of distribution of ImmuneScore, StromalScore, and ESTIMATEScore in the age classification, (D–F) in the FIGO stage, (G–I) in grades, (J–L) in T classifications, (M–O) in N classifications, (P–R) in M classifications.

The median TME score was used to divide patients into low and high TME groups, which can help us to identify a correlation between the TME and gene expression (Fig. 4). 1131 upregulated DEGs and 632 downregulated DEGs were acquired from ImmuneScore. 1299 upregulated DEGs and 51 downregulated DEGs were obtained from the StromalScore group. The intersection of the two groups, including upregulation and downregulation, was considered the final representation of the DEGs, which contained a total of 791 genes (749 upregulated genes and 42 downregulated genes). Moreover, these genes were used for the GO and KEGG enrichment analyses (Fig. 5). Based on the GO analysis, these DEGs are associated with immune functions, such as plasma membrane structure and leukocyte mediation (Fig. 5A–C). According to KEGG, DEGs are mainly involved in chemokine signaling pathways and cytokine–cytokine receptor interactions (Fig. 5D). It appears that immune factors play a pivotal role in cervical cancer TME.

Fig. 4.

Distribution of DEGs: Volcano plots, Venn diagrams, and Heatmap. Volcano plots showing the significance of DEGs between groups with high and low (A) ImmuneScore and (D) StromalScore. Venn diagrams showing the DEGs that were matched by the ImmuneScore and StromalScore assessments and were (B) upregulated and (E) downregulated. Heatmap for DEGs produced by contrasting the (C) ImmuneScore and (F) StromalScore groups with high and low scores. The gene name appears in the heatmap’s row, while the sample IDs that aren’t displayed in the plot appear in its column.

Fig. 5.

GO and KEGG enrichment analyses for 791 DEGs. (A–C) GO and (D) KEGG enrichment analyses, terms with p and q 0.05 were thought to be considerably enriched.

PPI networks, which were based on DEGs established in the STRING database and shown by Cytoscape software, were used to explore potential associations between these factors (Fig. 6). Subsequently, the top 25 key PPI genes (namely, PPI-key genes) were screened via the cytoHubba plugin in Cytoscape (Fig. 6A,B). 147 genes were significantly associated with the 791 DEGs using univariate Cox proportional hazard regression analysis (namely, COX-key genes). Thereafter, 11 genes existing in both groups of PPI- and COX-key genes were defined as being prognosis key genes and were visualized in the Venn diagram (Fig. 6C). Finally, we selected GRAP2, which was the best hazard ratio gene, for further analysis (Fig. 6D).

Fig. 6.

PPI network visualization with univariate Cox regression analyses. (A) An experimental evidence type was selected to build the PPI network in the STRING database. The size and color of the node represent the degree value, the color of the edge represents the total score value, and the edge represents the interaction between the target gene. The degree value that corresponds to the color from yellow to blue increases with the size of the node. (B) Bar-plot diagrams of the top 25 genes identified by Cytoscape’s cytoHubba plugin. (C) A Venn diagram showing where the important COX and PPI genes overlap. (D) The GRAP2 gene is highlighted in red on the Forest plot, which shows the substantial deferentially altered genes between the two cohorts.

The GRAP2 gene encodes GRB2-related adapter protein 2, which is a 37 kDa protein in humans and is also known as a GRB2-related adaptor downstream of Shc (GADS). GADS is involved in Ras signaling [17]. The latest research suggests that GADS is recruited to phosphorylated CD28 to activate T cells [18]. Increasing evidence has demonstrated that GRAP2 seems to serve as a protective factor in lung cancer [19]. Furthermore, in this study, we outlined the relationship between GRAP2 expression and clinical information in the TCGA-CESC database (Table 1) and observed that GRAP2 expression is closely associated with survival time and M stage in cervical cancer patients. Besides, the GRAP2 high-expression group possessed a relatively longer OS than the low-expression group (Fig. 7A). To avoid biased results due to a single database source, we selected 55 samples from GEO (GSE:52903) to verify this result. Similarly, we found that high GRAP2 expression predicted better OS (Fig. 7B). After conducting an analysis of GRAP2 combined with clinicopathological characteristics, we found that GRAP2 expression gradually decreased following the advanced FIGO stage and M classification (Fig. 7C–H). It is widely recognized that GRAP2 expression is associated with improved prognosis, including survival rates and metastasis rates, in cervical cancer.

Fig. 7.

Relationship between GRAP2 expression and clinical information. Low and high GRAP2 expression with (A) TCGA-CESC patients survival analyses (p 0.001, HR = 0.45, 95% CI = 0.27–0.73), and with (B) GEO (GSE:52903) patients (p = 0.044, HR = 0.41, 95% CI = 0.17–0.99). (C–H) Clinical information was associated with GRAP2 expression using either the Wilcoxon rank sum test or the Kruskal-Wallis rank sum test.

According to the abovementioned results, we believe that GRAP2expression may play a role in the malignant progression of cervical cancer. To further verify our theory, GSEA was used to explore the enrichment of GRAP2 expression in function and pathway. Genes within the GRAP2 high-expression tended to be associated with immune processes, such as cell activation involved in immune responses, cytokine-mediated signaling pathways, and mononuclear migration (Fig. 8A). Moreover, in the KEGG pathway analysis, the group mainly included immune-related pathways, including cytokine receptor interaction, chemokine signaling pathway, antigen processing and presentation, and T-cell receptor signaling pathway (Fig. 8C). Synchronously, the GRAP2 low-expression cohort was enriched in several different biological processes, including binding of sperm to the zona pellucida, establishment of mitochondrion localization, isoprenoid biosynthetic process, and negative regulation of nuclear division (Fig. 8B). Similarly, a wide variety of metabolic function pathways were enriched in low expression, which included the TCA cycle, glycerolipid metabolism, GPI anchor biosynthesis, and glyoxylate and dicarboxylate metabolism (Fig. 8D). In this figure, only the top ten items of pathways are shown in each section. All of the results illustrated that GRAP2 may serve as a promising indicator of immune and metabolic statuses in cervical cancer.

Fig. 8.

GSEA for GRAP2 expression. Up-regulated genes are situated on the left, which approaches the origin of the coordinates; in contrast, down-regulated genes are located on the right of the x-axis. Each line represents one specific gene set with a distinct color. The only gene sets that were deemed significant have NOM p 0.05 and FDR q 0.25. (A) Top 10 significant GOBP associated with high GRAP2 expression. (B) Top 10 significant pathways associated with low GRAP2 expression. (C) The enriched gene sets in KEGG with high GRAP2 expression. (D) The enriched gene sets in KEGG with the low GRAP2 expression.

The CIBERSORT algorithm was applied to further investigate the interplay between GRAP2 expression and TICs. First, Among the 22 immune cell types, 11 TICs had significant differences between the GRAP2 expression groups (Fig. 9A). Coincidentally, 11 TICs were strongly associated with GRAP2 expression from the correlation and difference analyses (Fig. 9B). Surprisingly, the intersection of the analyses showed a perfect coincidence of the two analyses (Fig. 9C). These TICs contain CD8 T cells, resting memory CD4 T cells, activated memory CD4 T cells, follicular helper T cells, M0 macrophages, M1 macrophages, activated dendritic cells, resting dendritic cells, resting mast cells, activated mast cells, eosinophils. According to these results, GRAP2 is essential to the immune TME of cervical cancer.

Fig. 9.

The proportions of TICs in groups with high and low GRAP2 expression and their correlation with GRAP2 expression. (A) By using the Wilcoxon rank sum test, 22 immune cell types in tumor tissues with low (green) or high (red) GRAP2 expression were compared. (B) An illustration of Pearson’s correlation between the number of TICs and GRAP2 expression, with blue lines denoting the best fit linear models. (C) An intersection between analyses in (A,B) shows 11 TICs shared by both.