HPV-Associated Gene Signatures in Bladder Cancer: A Comprehensive Prognostic Model and its Implications in Immunotherapy

Abstract

Background: Evidence increasingly indicates that HPV infection plays a pivotal role in the initiation and progression of bladder cancer (BC). Yet, determining the predictive value of HPV-associated genes in BC remains challenging. Methods: We identified differentially expressed HPV-associated genes of BC patients from the TCGA and GEO databases. We screened prognostic genes using COX and LASSO regression, subsequently establishing a risk prediction model. The model's precision and clinical relevance were gauged using Kaplan-Meier survival analyses and ROC curves. Functional enrichment, immune cell infiltration, and drug sensitivity analyses were performed across both high-risk and low-risk sets. PCR assays were utilized to measure the expression levels of genes. Results: We identified 13 HPV-associated genes for our risk model. Among these, FLRT2, HOXC5, LDLR, SCD, GRM7, DSC1, EMP1, and HMGA1 were identified as risk contributors, while LPA, SERPINA6, ZNF124, ETV7, and SCO2 were deemed protective. Cox regression analysis verified that our model provides an independent prediction of overall survival (OS) in bladder cancer (BC) patients. Gene Ontology (GO) analysis revealed predominant gene enrichment in wound healing, extracellular matrix composition, and collagen-rich extracellular matrices. KEGG pathway analysis highlighted primary enrichment areas, including focal adhesion, the PI3K-Akt signalling pathway, and ECM-receptor interaction. Risk scores were correlated with tumor microenvironment (TME) scores, immune cell infiltration, and sensitivities to both chemotherapy and immunotherapy. Conclusion: We have formulated a risk-assessment model pinpointing 13 central HPV-associated genes in BC. These genes present potential as prognostic indicators and therapeutic targets, emphasizing the intertwined relationship between HPV-induced BC progression and the immune landscape.

Keywords: Bladder cancer; Human papillomavirus (HPV); Prognostic model.

Figure 1

Heatmap showing HPV-related gene expression in BC samples versus normal samples. Horizontal coordinates represent samples, while vertical coordinates correspond to HPV-related genes represent.

Figure 2

(A, B) LASSO result diagram of TCGA set and (C, D) Risk curves and survival scatter plots of BC patients and (E) Heatmap of 13 HPV-related genes expression.

Figure 3

Evaluation of prognostic risk models for HPV-related genes. (A) Kaplan-Meier survival curves of two groups of BC patients. (B) Survival curves of PFS of two groups of BC patients. (C) Survival curves of BC patients during stage I-II. (D) Survival curves of BC patients during stage III-IV. (E) ROC curves of predictive efficacy to 1, 3, 5 years. (F) ROC curves of HPV-related risk model and other clinical prognostic factors.

Figure 4

(A) Univariate and (B) multivariate Cox regression analyses of the relationship between multiple clinical variables (including risk scores). (C) Kaplan-Meier survival curves for both groups of BC patients. (D) ROC curves to assess the predictive efficacy of the risk model. (E) Kaplan-Meier survival analysis of the H-TMB group and L-TMB group and (F) Kaplan-Meier survival analysis of the H-TMB high-risk group, H-TMB low-risk group, L-TMB high‑risk group and L-TMB low-risk group.

Figure 5

Heatmap of gene expression based on RNA-seq. (A) Between normal bladder cells and corresponding HPV-cells. (B) Between the wild-type bladder cancer (NC) strain and the bladder cancer cisplatin-resistant (CR) strain (C-F) Relative mRNA level adjusted to GAPDH of HMGA1, SCD, SCO2 and LDLR in SV-HUC-1, T24, J82, UM-UC3 and 5637 cell lines through RT-qPCR. ns ? **** ?? ***p < =0.001 **p <= 0.01 and *p <=0.05.

Figure 6

Relative mRNA level adjusted to GAPDH of DSC1, EMP1, ETV7, FLRT2, GRM7, HOXC5, LPA, SERPINA6 and ZNF124 in SV-HUC-1, T24, J82, UM-UC3 and 5637 cell lines through RT-qPCR. ***p < =0.001 **p <= 0.01 and *p <=0.05.

Figure 7

GO and KEGG analysis. (A, B, C) GO analysis on the biological processes (BP), cellular components (CC), and molecular functions (MF). (D, E) The KEGG pathway enrichment analysis.

Figure 8

(A-B) TME immunocyte infiltration characteristics and immune constituents in different risk groups. (A) Comparative analysis of stromal, immunity, and estimate scores between high-risk and low-risk groups. (B) Bar plot representing the differential proportions of 22 types of immune cell between high-risk and low-risk groups. (C-D) Immunocyte infiltration in different risk groups. (C) Boxplot for the immune-associated functions. (D) Boxplot of the immune cell abundance. (*p < 0.05; ** p < 0.01; *** p < 0.001) E-H Univariate and multivariate Cox regression analysis of the correlation between OS and various clinical variables (E). Kaplan-Meier survival analysis (F) and correlation between the risk score and clinical response (complete response [CR]/partial response [PR] and stable disease [SD]/progressive disease [PD]) (G) and ROC curve (H) in the IMvigor210 cohort. ***p < =0.001 **p <= 0.01 and *p <=0.05.