We downloaded three datasets of gene expression profiles in cervical cancer, including TCGA-CESC (5), GSE29570 (11), and GSE39001 (12). The TCGA-CESC dataset was from the genomic data commons (GDC) data portal (https://portal.gdc.cancer.gov/), and GSE30784 and GSE39366 were from the NCBI gene expression omnibus (https://www.ncbi.nlm.nih.gov/geo/). From GDC, we also downloaded the protein expression, somatic mutation, and somatic copy number alteration (CNA) profiles and clinical data in TCGA-CESC. We took TCGA-CESC as the discovery set and performed main analyses in this dataset. The other two datasets were validation sets. These datasets contained only data of human mRNA expression, but not data of the expression of infecting HPV variants. The data on the mRNA expression of HPV alpha-7 and alpha-9 clades were not available in any of these datasets. A summary of these datasets is shown in Table 1 .

We first identified 50 genes having the largest expression variations across the HPV+ cervical cancers in TCGA-CESC. Based on the expression profiles of the 50 genes, we performed the hierarchical clustering of the HPV+ cervical cancers in the three datasets, respectively.

We quantified the enrichment level of an immune signature or phenotypic feature in a tumor sample by the single-sample gene-set enrichment analysis (ssGSEA) of its marker gene set (13). The ssGSEA calculates the enrichment score of a gene set in a sample based on its expression profiles. The ratios of immunostimulatory/immunosuppressive signatures were the base-2 log-transformed values of the geometric mean expression levels of all marker genes of immunostimulatory signatures divided by those of immunosuppressive signatures. The marker gene sets of immune signatures or phenotypic features are shown in Supplementary Table S1 .

We first identified the genes differentially expressed between two classes of samples using the Student’s t test with a threshold of adjusted P value < 0.05 and fold change of their mean expression levels > 2. By inputting the differentially expressed genes into the GSEA web tool (14), we identified KEGG (15) pathways highly enriched in one class versus another class using a threshold of adjusted P value < 0.05. We used the weighted gene co-expression network analysis (WGCNA) (16) to identify gene modules significantly enriched in subtypes. Based on the expression correlations between the gene modules’ hub genes, we identified GO terms having significant correlations with specific trait by WGCNA.

We used the ESTIMATE algorithm (17) to calculate tumor immune score, stromal score, and tumor purity based on gene expression profiles. The immune score, stromal score, and tumor purity indicate the level of tumor-infiltrating lymphocytes, proportion of stromal component, and proportion of tumor cells in a bulk tumor sample. The stemness score was calculated by the ssGSEA of its marker gene set (18) and defined the level of tumor stem cell-like phenotypic feature. We evaluated the ITH level by the DEPTH algorithm (19), which quantifies ITH based on transcriptomic alterations in the tumor.

We used Kaplan–Meier curves to display the survival time difference between cervical cancer subtypes, whose significance was evaluated by the log-rank test. The function “survfit” in the R package “survival” was utilized to perform this analysis. In addition, we performed multivariate survival analysis to investigate the correlation between cervical cancer subtypes and survival prognosis after adjusting for confounding variables, including immune score, stemness score, and ITH score. All these variables are continuous variables. The Cox proportional hazards model was utilized to implement the multivariate survival analysis with the function “coxph” in the R package “survival”.

The TMB was defined as the total number of somatic mutations in the tumor. We used GISTIC2 (20) to calculate G-scores for cervical cancer subtypes. The G-score represents the amplitude of the CNA and the frequency of its occurrence across a group of samples (20).

We utilized the Random Forest (RF) algorithm (21) to predict cervical cancer subtypes. In the RF, the number of trees was set to 500, and the features included the 50 genes with the largest expression variations across the HPV+ cervical cancers in TCGA-CESC. We reported sensitivity, specificity, and area under the receiver operating characteristic curve (AUC) to evaluate the prediction performance. The R package “randomForest” was utilized to perform the RF algorithm. The R package “PreHPVcc” for predicting HPV+ cervical cancer subtypes is available in the GitHub repository (https://github.com/WangX-Lab/PreHPVcc).

In class comparisons, we used the Mann–Whitney U test or Kruskal–Wallis (K–W) test for not normally distributed data (Shapiro test, P < 0.05) and Student’s t test or one-way analysis of variance (ANOVA) test for normally distributed data. We utilized the chi-square test to analyze contingency tables. We used the Benjamini-Hochberg method (22) to adjust for P values in multiple tests. We performed all statistical analyses in the R programming environment (version 4.0.2).

We hierarchically clustered HPV+ cervical cancers based on the expression profiles of 50 genes which had the largest expression variations across the HPV+ cervical cancers. We obtained two clear clusters, termed HPV+G1 and HPV+G2 ( Figure 1A ). We confirmed that this classification was reproducible in two other datasets (GSE29570 and GSE39001) ( Figure 1A ). We found that HPV+G1 had a significantly higher disease-free survival (DFS) rate than HPV+G2 (log-rank test, P = 0.01) ( Figure 1B ). Interestingly, HPV+G1 showed significantly higher enrichment levels of various immune signatures than HPV+G2, including CD8+ T cells, B cells, M1 macrophages, cytolytic activity, IFN response, CD4+ regulatory T cells, pro-inflammatory cytokines, T cell exhaustion, MDSCs, PD-L1 expression, and anti-inflammatory cytokines (one-tailed Mann–Whitney U test, P < 0.01) ( Figure 1C ). Moreover, the ratios of immunostimulatory/immunosuppressive signatures (M1/M2 macrophages and pro/anti-inflammatory cytokines) were higher in HPV+G1 than in HPV+G2 ( Figure 1C ). We further used the ESTIMATE algorithm (17) to calculate the immune score representing the tumor immune infiltration level. As expected, immune scores were significantly higher in HPV+G1 than in HPV+G2 (P < 0.001) ( Figure 1D ), while tumor purity was significantly lower in HPV+G1 than in HPV+G2 (P < 0.001) ( Figure 1E ). These results suggest that HPV+G1 has a more active tumor immune microenvironment than HPV+G2. In addition, we found that stromal scores were significantly higher in HPV+G1 than in HPV+G2 (P = 0.002) ( Figure 1F ). HPV+G1 had significantly lower stemness scores than HPV+G2 (P = 0.005) ( Figure 1G ). Moreover, HPV+G1 had significantly lower ITH scores than HPV+G2 (P = 0.005) ( Figure 1H ).

Because tumor immune signatures, stemness, and ITH are associated with clinical outcomes in cancer (18, 23, 24) and had significantly different scores between the HPV+ cervical cancer subtypes, the survival difference between the subtypes could be impacted by these confounding variables. To explore the possibility, we performed multivariate (immune score, stemness score, ITH score, and subtype) survival analysis with the multivariate Cox proportional hazards model. The result showed that the subtype HPV+G2 remained a significant risk factor (P = 0.007; hazard ratio (HR) = 3.76 and its 95% confidence interval (CI) (3):) ( Figure 1I ). It suggests that the DFS difference between both subtypes is independent of these confounding variables.

Genomic instability often leads to increased TMB and CNAs (25). We found that HPV+G1 had higher TMB than HPV+G2 (one-tailed Mann–Whitney U test, P = 0.057) ( Figure 2A ). Homologous recombination deficiency (HRD) may lead to aneuploidy, namely CNAs (25). We obtained the HRD scores for the TCGA cervical cancers from the publication by Knijnenburg et al., which were the combined scores of HRD loss of heterozygosity, large-scale state transitions, and the number of telomeric allelic imbalances (25). We found that HPV+G1 had significantly higher HRD scores than HPV+G2 (one-tailed Mann–Whitney U test, P = 0.015) ( Figure 2B ). We found that the G-scores of copy number amplifications and deletions were significantly higher in HPV+G1 than in HPV+G2 cervical cancers ( Figure 2C ). Taken together, these results indicated that HPV+G1 was more genomically instable than HPV+G2.

We compared global methylation levels (26) between both subtypes and found that HPV+G1 had significantly lower global methylation levels than HPV+G2 (one-tailed Mann–Whitney U test, P = 0.011) ( Figure 2D ). This result conforms with that low methylation levels is associated with increased TMB and CNAs in cancer (26). Strikingly, we found that 5367 genes showed significantly lower methylation levels in HPV+G1 than in HPV+G2 (FDR < 0.05), while there was no any gene showing significantly higher methylation levels in HPV+G1 than in HPV+G2. These results indicate a significant difference in epigenomic profiles between both subtypes.

We compared gene expression profiles between HPV+G1 and HPV+G2 and identified significantly upregulated genes in both subtypes. Based on these upregulated genes, we identified KEGG pathways highly enriched in HPV+G1 and HPV+G2, respectively, using the GSEA web tool (14). Many of the pathways especially enriched in HPV+G1 were involved in immune signatures, including cytokine-cytokine receptor interaction, cell adhesion molecules, complement and coagulation cascades, chemokine signaling, Toll-like receptor signaling, Fc gamma R-mediated phagocytosis, leukocyte transendothelial migration, intestinal immune network for IgA production, natural killer cell-mediated cytotoxicity, NOD-like receptor signaling, and Jak-STAT signaling ( Figure 3A ). It confirms the more active tumor immune microenvironment in HPV+G1 versus HPV+G2.

WGCNA (16) identified nine gene modules significantly differentiating cervical cancers by subtype ( Figure 3B ). The representative GO terms for the gene modules upregulated in HPV+G1 while downregulated in HPV+G2 included immune response (in brown) and epidermal and epithelial cell differentiation (in black). In contrast, the representative GO terms for the gene modules upregulated in HPV+G2 while downregulated in HPV+G1 included intracellular non-membrane-bounded organelle (in green) and microvillus (in pink). Again, these results confirm that HPV+G1 has a more active tumor immune microenvironment versus HPV+G2.

We compared the expression levels of 192 proteins between HPV+G1 and HPV+G2 and identified significantly upregulated proteins in both subtypes (two-tailed Student’s t test, P < 0.05). We found 22 proteins having significantly higher expression levels in HPV+G1, including NDRG1_pT346, Notch1, EGFR, Annexin-1, CD49b, PI3K-p85, EGFR_pY1068, Caveolin-1, YB-1_pS102, Src_pY416, Bad_pS112, MEK1, Myosin-IIa_pS1943, PAI-1, YAP_pS127, YAP, Bcl-2, Fibronectin, MYH11, GSK3_pS9, Syk, and VHL ( Figure 4 ). Among them, Annexin A1, also known as lipocortin I, plays a role in the regulation of innate and adaptive immunity (27). CD49b (cluster of differentiation 49b) is an integrin alpha subunit expressed on various cell types, including T cells and NK cells (28). YAP is involved in the Hippo signaling pathway, which plays a role in tumor immune regulation (29). Bcl-2 is a positive regulator of apoptosis, which is associated with anti-tumor immunity (30). Collectively, these results again confirmed the more active tumor immune microenvironment in HPV+G1 versus HPV+G2.

In contrast, 25 proteins showed significantly higher expression levels in HPV+G2 than in HPV+G1 ( Figure 4 ). These proteins included TIGAR, XRCC1, AMPK-α, Claudin-7, Bcl-xL, PCNA, p27, SCD1, TAZ, INPP4B, Ku80, ASNS, AMPK_pT172, MSH2, Acetyl-a-Tubulin-Lys40, 53BP1, Rab25, MIG-6, MSH6, mTOR, Chk2, ER-α, c-Kit, PKC-delta_pS664, and Transglutaminase. Among them, many proteins are involved in DNA repair, including XRCC1, PCNA, Ku80, MSH2, and MSH6, confirming that HPV+G2 was more genomically stable than HPV+G2.

Although HPV+G2 had lower enrichment levels of immune signatures than HPV+G1, it showed significantly higher enrichment levels of various immune signatures than HPV- tumors, including

NK cells, M1 macrophages, IFN response, CD4+ regulatory T cells, and M2 macrophages ( Supplementary Figure S1A ). Moreover, the ratios of immunostimulatory/immunosuppressive signatures (CD8+ T cells/MDSCs) were significantly higher in HPV+G2 than in HPV- cervical cancers ( Supplementary Figure S1B ). Most of the human leukocyte antigen (HLA) genes, which encode major histocompatibility complex (MHC) proteins and play essential roles in the regulation of the immune system, displayed significantly higher expression levels in HPV+G2 than in HPV- cervical cancers (two-tailed Student’s t test, FDR < 0.02; fold change (FC) > 1.5) ( Supplementary Figure S1C ). Immune scores were significantly higher in HPV+G2 than in HPV- cervical cancers (P = 0.04) ( Supplementary Figure S1D ). These results support that HPV infection results in a more active tumor immune microenvironment in cervical cancer.

Both ITH and stemness scores followed the pattern: HPV+G1 < HPV+G2 < HPV-, while global methylation levels were higher in HPV+G2 than in both HPV+G1 and HPV- ( Supplementary Figure S1E ). In addition, although TMB and HRD scores were higher in HPV+G1 than in HPV+G2, they were not significantly different between HPV+G2 and HPV- (P > 0.2). Furthermore, we did not observe significantly different OS or DFS rate between the HPV+ subtypes and HPV- (log-rank test, P > 0.2).

We used TCGA-CESC as the training set and the other two datasets as test sets. The 10-fold cross-validation (CV) sensitivity, specificity, and AUC in TCGA-CESC were 99.1%, 95.0%, and 100.0%, respectively. The prediction sensitivity, specificity, and AUC in GSE29570 were 100%, 80.0%, and 90.0%, respectively, and in GSE39001 were 92.0%, 100%, and 96.0%, respectively ( Figure 5 ). These results suggest that our subtyping method for HPV+ cervical cancers is highly reproducible and predictable. In the prediction model, we found that the 10 features (genes) with the highest importance weights included DSG3, DSC3, CLCA2, LASS3, CALML3, SERPINB13, IVL, PROM1, AGR3, and TMPRSS11D ( Table 2 ).

We investigated the relationship between our subtyping method and other cervical cancer subtyping methods (5). We found that squamous cell carcinomas constituted 98% of HPV+G1 tumors versus 34% of HPV+G2 tumors (chi-square test, P < 0.001) ( Figure 6 ). In contrast, adenocarcinomas constituted 63% of HPV+G2 tumors versus 2% of HPV+G1 tumors. The TCGA network classified cervical cancers into three subtypes by unsupervised clustering of variable DNA-methylation probes (5). The three subtypes included: CpG island hypermethylated (CIMP-high), CIMP-intermediate, and CIMP-low. We found that 4% of HPV+G1 tumors were CIMP-high, compared to 53% of HPV+G2 tumors being CIMP-high, and that 46% of HPV+G1 tumors were CIMP-low versus 29% of HPV+G2 tumors being CIMP-low (P < 0.001) ( Figure 6 ). It is consistent with the lower overall DNA methylation level in HPV+G1 relative to HPV+G2. In addition, we found 20% of HPV+G1 tumors being Clade A7 versus 50% of HPV+G2 tumors and 78% of HPV+G1 tumors being Clade A9 versus 50% of HPV+G2 tumors (P = 0.002) ( Figure 6 ). The different distribution of HPV clades between HPV+G1 and HPV+G2 indicates a better prognosis in HPV+G1 versus HPV+G2 since HPV clade A7 cervical cancers are more aggressive than clade A9 cancer (31). RPPA-based clustering identified three clusters: hormone, EMT, and PI3K-AKT (5). We found that 34% of HPV+G1 tumors were in the EMT cluster versus 8% of HPV+G2 tumors (P = 0.002) ( Figure 6 ). It is justified since EMT represents a stromal signature and HPV+G1 has higher stromal scores than HPV+G2. In addition, 32% of HPV+G1 tumors were in the PI3K-AKT cluster versus 17% of HPV+G2 tumors, and 34% of HPV+G1 tumors were in the hormone cluster versus 75% of HPV+G2 tumors. It is consistent with previous results that PI3K-p85 was more abundant in HPV+G1 while ER-α was more abundant in HPV+G2. Furthermore, 37% of HPV+G1 tumors were HPV-16 positive versus 32% of HPV+G2 tumors, and 5% of HPV+G1 tumors were HPV-18 positive versus 25% of HPV+G2 tumors. It indicates that a significantly higher proportion of HPV+G2 tumors are HPV-18 positive compared to HPV+G1 tumors (P < 0.001). This result supports that HPV-18 infection is an adverse prognostic factor in cervical cancer (3), while HPV-16 infection is a positive prognostic factor in cervical cancer (32).