Identification of Tumor Microenvironment and DNA Methylation-Related Prognostic Signature for Predicting Clinical Outcomes and Therapeutic Responses in Cervical Cancer

Abstract

Background: Tumor microenvironment (TME) has been reported to have a strong association with tumor progression and therapeutic outcome, and epigenetic modifications such as DNA methylation can affect TMB and play an indispensable role in tumorigenesis. However, the potential mechanisms of TME and DNA methylation remain unclear in cervical cancer (CC). Methods: The immune and stromal scores of TME were generated by the ESTIMATE algorithm for CC patients in The Cancer Genome Atlas (TCGA) database. The TME and DNA methylation-related genes were identified by the integrative analysis of DNA promoter methylation and gene expression. The least absolute shrinkage and selection operator (LASSO) Cox regression was performed 1,000 times to further identify a nine-gene TME and DNA methylation-related prognostic signature. The signature was further validated in Gene Expression Omnibus (GEO) dataset. Then, the identified signature was integrated with the Federation International of Gynecology and Obstetrics (FIGO) stage to establish a composite prognostic nomogram. Results: CC patients with high immunity levels have better survival than those with low immunity levels. Both in the training and validation datasets, the risk score of the signature was an independent prognosis factor. The composite nomogram showed higher accuracy of prognosis and greater net benefits than the FIGO stage and the signature. The high-risk group had a significantly higher fraction of genome altered than the low-risk group. Eleven genes were significantly different in mutation frequencies between the high- and low-risk groups. Interestingly, patients with mutant TTN had better overall survival (OS) than those with wild type. Patients in the low-risk group had significantly higher tumor mutational burden (TMB) than those in the high-risk group. Taken together, the results of TMB, immunophenoscore (IPS), and tumor immune dysfunction and exclusion (TIDE) score suggested that patients in the low-risk group may have greater immunotherapy benefits. Finally, four drugs (panobinostat, lenvatinib, everolimus, and temsirolimus) were found to have potential therapeutic implications for patients with a high-risk score. Conclusions: Our findings highlight that the TME and DNA methylation-related prognostic signature can accurately predict the prognosis of CC and may be important for stratified management of patients and precision targeted therapy.

Keywords: DNA methylation; cervical cancer; drug response; immunotherapy response; prognostic model; tumor microenvironment.

FIGURE 1

Identification of the TME and DNA methylation-related prognostic signature. (A) Scatter plot of promoter mean methylation difference and gene expression levels change. hyper-up, hypermethylated-upregulated; hyper-down, hypermethylated-downregulated; hypo-up, hypomethylated-upregulated; hypo-down, hypomethylated-downregulated. (B,C) Venn diagrams showing the intersection between DEGs and hypermethylated genes (top) and between DEGs and hypomethylated genes (bottom). (D) The C-index of different genes combinations in the signature. (E) The nine genes included in the signature. Corresponding coefficients are depicted by horizontal bars.

FIGURE 2

Validation of the prognostic value of the risk score. (A) Difference analysis of the distribution of risk scores in different FIGO stages, TNM stages, and histological types. (B) Kaplan–Meier curves for differential detection of patients in the TCGA cohort by the log-rank test. (C) ROC curves of risk scores used to predict 1-year, 3-year, and 5-year survival in the TCGA cohort. (D) Time-dependent ROC curves of the risk score in the TCGA cohorts. (E) Kaplan–Meier curves for differential detection of patients in the GSE44001 cohort by the log-rank test. (F) ROC curves of risk scores used to predict 1-year, 3-year, and 5-year survival in the GSE44001 cohort. (G) Time-dependent ROC curves of the risk score in the GSE44001 cohorts. *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 3

Forest plot of the univariate and multivariate Cox regression analysis in TCGA and GSE44001 cohorts.

FIGURE 4

Construction of a nomogram model. (A) Nomogram constructed in conjunction with the risk score and FIGO stage for the TCGA cohort. (B) Calibration plot of the nomogram. (C) C-index curves of the FIGO stage, risk score, and nomogram. (D) Decision curve analysis for evaluating the net benefits of FIGO stage, risk score, and nomogram.

FIGURE 5

The immune signature between the high- and low-risk groups in the TCGA cohort. (A–C) Association between immune score, stromal score, tumor purity, and risk score and their distribution in different risk groups. (D,E) Differential analysis of gene expression levels of HLA family genes and immune checkpoints in different risk groups. (F) Correlation analysis for the risk score and the gene expression levels of HLA family genes and immune checkpoints. (G) The heatmap showing the immune and stromal cell infiltration levels and differences in distribution between different risk groups. *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 6

Function analysis of genes correlated with the risk score. (A) GSEA enrichment plots showing enriched gene sets against to hallmark dataset in high- and low-risk groups. NES, normalized enrichment score. (B) A dendrogram of the top 5,000 genes with the most variation clustered based on the topological overlap together. (C) The heatmap showing the association between gene modules and the signature risk score. (D) GSEA annotated by KEGG gene sets for the brown module genes.

FIGURE 7

Identification differences of the genetic variation and pathway activation between high- and low-risk groups. (A) Tumor mutation burdens were compared among distinct risk groups. (B) Forest plot of genes with differences in mutation frequencies between the low- and high-risk groups. (C) Waterfall plot of the 11 mutant genes with significant frequency differences between low- and high-risk groups. (D) Interaction of differentially mutated genes. (E) Kaplan–Meier curve showing that patients with mutant TTN have a better OS than those with wild type. (F) Differential analysis of GSVA scores among distinct risk groups. (G) Copy number alteration gains (red) and losses (blue) between the low- and high-risk groups. (H) Differential analysis of altered, lost, and gained genome fractions (%) between the low-risk and high-risk groups. *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 8

Identification of potential agents and prediction of immunotherapeutic effect. (A,B) Differential drug response analysis of the selected agents for CC patients between the higher and lower risk score groups based on the CTRP dataset and Spearman’s correlation analysis of CTRP-derived agents and risk score. (C,D) Differential drug response analysis of the selected agents for CC patients between the higher and lower risk score groups based on the PRISM dataset and Spearman’s correlation analysis of PRISM-derived agents and risk score. (E,F) The TIDE score and IPS were compared between the high- and low-risk groups. *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.