Integrating bioinformatics to explore HPV-31 and HPV-52 E6/E7 proteins: from structural analysis to antigenic epitope prediction

TL;DR

This study uses bioinformatics to analyze the E6 and E7 proteins of HPV-31 and HPV-52, identifying key epitopes that could help in developing vaccines against cervical cancer.

Contribution

The study provides a detailed, subtype-specific analysis of immunogenic epitopes in HPV-31 and HPV-52 E6/E7 proteins.

Findings

HPV-31 and HPV-52 E6/E7 proteins have distinct B-cell and T-cell epitopes identified through bioinformatics analysis.

The instability indices of the E6 and E7 proteins suggest structural instability, which may influence their function and antigenicity.

The findings offer a foundation for experimental validation and vaccine development targeting HPV-31 and HPV-52.

Abstract

Cervical cancer is the most common malignant neoplasm of the female reproductive tract. Infection with human papillomavirus (HPV) has been strongly associated with cervical cancer. Previous bioinformatics studies have examined the E6 and E7 proteins of high-risk HPV types; however, subtype-specific analyses for HPV-31 and HPV-52 remain limited. Understanding the structure and properties of the E6 and E7 proteins of HPV-31 and HPV-52 is crucial to elucidating their functions and advancing vaccine development.

A bioinformatics approach was employed to predict the physicochemical properties, hydrophilicity, protein structure, glycosylation sites, phosphorylation sites, terminal positions, signal peptide cleavage sites, transmembrane regions, homology, and dominant epitopes of the E6 and E7 proteins of HPV-31 and HPV-52.

For HPV-31 E6, an instability index (II) of 43.93 indicated that the protein is unstable; potential B-cell epitopes were identified at residues 55–61 (RDDTPYG), 112–116 (PEEKQ), and 125–131 (FHNIGGR), while T-cell epitopes were predicted at residues 45–53 (FAFTDLTIV) and 72–80 (KVSEFRWYR). HPV-52 E6 exhibited an instability index (II) of 55.57, with B-cell epitopes at residues 110–119 (LCPEEKERHV) and 129–141 (MGRWTGRCSECWR), and T-cell epitopes at residues 45–53 (FLFTDLRIV) and 82–87 (SLYGKT). HPV-31 E7, with an instability index (II) of 51.05, exhibited B-cell epitopes at residues 8–17 (QDYYLDLQP), 16–20 (QPEAT), 29–41 (PDSSDEEDVIDEP), and 42–48 (AGQAKPDT), and T-cell epitopes at residues 7–15 (TLQDYVLDL) and 82–90 (LLMGSFGIV). HPV-52 E7, with an instability index (II) of 49.15, exhibited B-cell epitopes at residues 11–19 (YILDLQPET), 23–27 (HCYEQ), 29–38 (GDSSDEEDTD), and 36–48 (DTDGVDRPDGQAE), and T-cell epitopes at residues 53–59 (NYYIVTY) and 84–90 (MLLGTLQ).

In summary, the E6 and E7 proteins of HPV-31 and HPV-52 contain dominant epitopes for both T cells and B cells. These findings delineate subtype-specific immunogenic regions and establish a foundation for experimental validation and vaccine design.

Figures

Phosphorylation sites: (A) HPV31 E6 (B) HPV31 E7 (C) HPV52 E6 (D) HPV52 E7.

Phosphorylation sites: (A) HPV31 E6 (B) HPV31 E7 (C) HPV52 E6 (D) HPV52 E7.

TMHMM analyzed the transmembrane domain of the proteins. (A) HPV31 E6 (B) HPV31 E7 (C) HPV52 E6 (D) HPV52 E7.

Secondary structure prediction: (A) HPV31 E6 oncoprotein; (B) HPV31 E7 oncoprotein; (C) HPV52 E6 oncoprotein; (D) HPV52 E7 oncoprotein.

Tertiary structure prediction. (A) HPV31 E6 protein; (B) HPV31 E7 protein; (C) HPV52 E6 protein; (D) HPV52 E7 protein.

Homology and molecular evolution analysis. (A) Homology analysis of E6 proteins of HPV; (B) The molecular evolutionary tree of E6 proteins of HPV; (C) Homology analysis of E7 proteins of HPV; (D) The molecular evolutionary tree of E7 proteins of HPV.

Tables

Table 1. Methods summary table.

Step	Tool/Method used	Purpose	Key parameters
Amino Acid Sequence	NCBI Database	Retrieve full-length protein sequences	HPV31 E6/E7, HPV52 E6/E7
Physicochemical Parameters	ProtParam, ProtScale	Calculate molecular weight, pI, hydrophobicity, etc.	ProtParam: instability index, GRAVY; ProtScale: hydrophobicity (Kyte–Doolittle), hydrophilicity (Hopp–Woods)
PTM Site Prediction	NetPhos 3.1, MotifScan, NetNGlyc	Predict phosphorylation, kinase motifs, glycosylation	Threshold: 0.5
Signal Peptide Prediction	SignalP 4.1	Predict signal peptide cleavage	D-score ≥ 0.45
Transmembrane Helix	TMHMM 2.0	Predict transmembrane regions	Probability ≥ 0.5
Secondary Structure	SOPMA v3.0	Predict secondary structure (α-helix, β-sheet, etc.)	Threshold: 8% difference, window width = 17
Tertiary Structure	Phyre2 v2.0	Homology modeling	Templates: ≥ 90% coverage, ≥ 99% confidence
Sequence Homology and Phylogenetic Analysis	Clustal X 2.0, MEGA 7.0.20	Align sequences and infer phylogenetic tree	Gap open penalty = 10, gap extension = 0.1, Bootstrap=1,000
B-cell Epitope Prediction	ABCpred v2.0, BepiPred 1.0, BCPREDS 1.0, SVMTrip v1.0	Predict linear epitopes for B-cells	ABCpred v2.0: threshold 0.51; peptide length = 16, BepiPred 1.0: threshold 0.35; window = 20, BCPREDS 1.0: epitope length = 20; specificity = 75%,SVMTrip v1.0: threshold 0.51; peptide length = 20
T-cell Epitope Prediction	SYFPEITHI v1.0, IEDB MHC II, IEDB MHC I	Predict CD4+ and CD8+ T-cell epitopes	CD4+: SYFPEITHI score ≥ 20; CD8+: NetMHCpan percentile rank ≤ 1

Table 2. Summarizes physicochemical parameters for all four proteins.

Protein	AA Length	Molecular Mass (Da)	Theoretical pI	Basic (K,R)	Acidic (D,E)	Instability Index (II)	GRAVY	Classification (II > 40 = unstable; GRAVY < 0 = hydrophilic)
HPV-31 E6	149	17,767.61	9.13	27	18	43.93	–0.567	Unstable; Hydrophilic ( Figure 1A )
HPV-31 E7	98	10,944.27	3.90	5	16	51.05	–0.235	Unstable; Hydrophilic ( Figure 1B )
HPV-52 E6	148	17,925.85	8.96	26	19	55.57	–0.599	Unstable; Hydrophilic ( Figure 1C )
HPV-52 E7	99	11,032.24	4.33	7	17	49.15	–0.459	Unstable; Hydrophilic ( Figure 1D )

Table 3. Summary of predicted PTM sites and membrane localization (NetPhos 3.1; MotifScan v2022; NetNGlyc 1.0; SignalP 4.1; TMHMM 2.0).

Protein	Phosphorylation (NetPhos > 0.5)	CK2 (CK2 motif > 0.5)	PKC (PKC motif > 0.5)	TK (TK motif > 0.5)	N-Glycosylation (NetNGlyc > 0.5)	Signal Peptide (SignalP D > 0.45)	Transmembrane (TMHMM > 0.5) (Figure 3)
HPV-31 E6	S16, S17, S71, S74, S82; T38, T40, T58, T64, T86, T133, T145, T147; Y60 ( Figure 2A )	17–20, 38–42, 86–89	92–94, 133–135	72–79	None	None	None
HPV-52 E6	S22, S71, S74, S82, S97; T11, T48, T108, T133, T146; Y60 ( Figure 2C )	11–14, 22–25, 87–90	100–102, 133–135	72–79	None	None	None
HPV-31 E7	S31, S32, S40, S50, S86; T5, T20, T64, T72; Y52 ( Figure 2B )	7–10, 31–34, 72–75	64–66, 95–97	None	None	None	None
HPV-52 E7	S31, S32; T7, T19, T20, T37, T58, T66, T76; Y11 ( Figure 2D )	7–10, 31–34, 74–77	7–9, 66–68	None	None	None	None

Table 4. Summarizes secondary structure content.

Protein	α-Helix (%)	β-Sheet (%)	β-Turn (%)	Random Coil (%)	Interpretation
HPV-31 E6	49.66	14.56	4.43	35.44	Relatively compact, fewer surface coils
HPV-31 E7	25.51	22.45	0.00	52.04	More random coils, implies greater surface exposure
HPV-52 E6	54.05	10.81	1.35	33.78	Compact with predominant α-helices
HPV-52 E7	27.27	21.21	0.00	51.52	Loose structure with significant coils

Table 5. HPV-31/52 E6/E7 B-Cell epitope candidates (ABCpred; BepiPred; BCPREDS; SVMTrip).

Protein	Algorithm Combination	Residues	Sequence	Structural Context	Tool Version/Threshold
HPV-31 E6	ABCpred & BepiPred & BCPREDS & SVMTrip	55–61	RDDTPYG	Random coil adjacent to LxxLL binding pocket	ABCpred v2.0 (threshold = 0.51)BepiPred 1.0 (threshold = 0.35)BCPREDS 1.0 (length = 20; specificity = 75%)SVMTrip v1.0 (length = 20; threshold = 0.51)
	BepiPred & BCPREDS	112–116	PEEKQ	β-turn at surface	Same as above
	ABCpred & SVMTrip	125–131	FHNIGGR	C-terminal random coil near functional region	Same as above
HPV-52 E6	ABCpred & BepiPred	110–119	LCPEEKERHV	C-terminal β-turn near Zn-finger	Same as above
	BepiPred & SVMTrip	129–141	MGRWTGRCSECWR	Random coil loop; structurally exposed	Same as above
HPV-31 E7	ABCpred & BepiPred	8–17	QDYYLDLQP	N-terminal random coil, high hydrophilicity	Same as above
	BepiPred & BCPREDS	16–20	QPEAT	Small β-turn in central region	Same as above
	BCPREDS & SVMTrip	29–41	PDSSDEEDVIDEP	Long random coil loop with high immunogenicity	Same as above
	SVMTrip & ABCpred	42–48	AGQAKPDT	C-terminal loop region accessible to antibodies	Same as above
HPV-52 E7	ABCpred & BepiPred	11–19	YILDLQPET	N-terminal random coil loop	Same as above
	BepiPred & BCPREDS	23–27	HCYEQ	Small β-turn	Same as above
	BCPREDS & SVMTrip	29–38	GDSSDEEDTD	Central random coil loop	Same as above
	SVMTrip & ABCpred	36–48	DTDGVDRPDGQAE	C-terminal loop region	Same as above

Table 6. HPV-31/52 E6/E7 T-Cell Epitope Candidates (SYFPEITHI; IEDB).

Protein	Type	HLA Allele	Residues	Sequence	Affinity Metric (IEDB percentile)	SYFPEITHI Score	Structural Context	Tool Version/Threshold
HPV-31 E6	CD4+	DRB1*1501	45–53	FAFTDLTIV	3.10	25	Zn-finger region; likely helper epitope	SYFPEITHI v1.0 (score ≥ 20)IEDB MHC-II (percentile ≤ 10)IEDB MHC-I (NetMHCpan 4.1; percentile ≤ 1)
	CD4+	DRB1*1501	72–80	KVSEFRWYR	4.50	22	β-turn at surface; T_H potentiation	Same as above
	CD8+	A*1101	82–90	SVYGTTLEK	0.01	—	Conserved helix; cross-subtype CTL potential	Same as above
	CD8+	A*0201	45–53	FAFTDLTIV	0.93	—	Overlaps with CD4+ 45–53; candidate for poly-epitope design	Same as above
HPV-52 E6	CD4+	DRB1*1501	45–53	FLFTDLRIV	3.70	24	Conserved block; cross-protection candidate	Same as above
	CD4+	DRB1*1501	82–87	SLYGKT	4.00	23	Loop region; potential helper epitope	Same as above
	CD8+	A*1101	86–94	KTLEERVRK	0.01	—	Zn-finger adjacency; strong CTL candidate	Same as above
	CD8+	A*0201	18–26	VLEESVHEI	0.03	—	N-terminal helix; antigen-presenting potential	Same as above
HPV-31 E7	CD4+	DRB1*1501	7–15	TLQDYVLDL	7.10	21	N-terminal random coil; T_H epitope candidate	Same as above
	CD4+	DRB1*1501	82–90	LLMGSFGIV	8.00	20	C-terminal coil; possible cross-reactive	Same as above
	CD8+	A*0201	7–15	TLQDYVLDL	0.09	—	Overlaps CD4+ 7–15; poly-epitope design potential	Same as above
	CD8+	A*1101	37–46	VIDSPAGQAK	0.33	—	Central coil loop; strong CTL candidate	Same as above
HPV-52 E7	CD4+	DRB1*1501	11–19	YILDLQPET	9.00	19	Loop region; intermediate T_H immunogenicity	Same as above
	CD4+	DRB1*1501	84–90	MLLGTLQ	13.00	18	C-terminal coil; modest helper response	Same as above
	CD8+	A*0201	84–92	MLLGTLQVV	0.08	—	C-terminal coil; strong CTL candidate	Same as above
	CD8+	A*1101	51–59	TSNYYIVTY	0.74	—	Central coil; potential CTL memory locater	Same as above

Full text

Introduction

1

Human papillomavirus (HPV) is among the most prevalent sexually transmitted viruses worldwide, and infection with HPV has been strongly associated with the development of various cancers, particularly cervical cancer (1). Since the landmark identification of HPV’s role in cervical carcinogenesis in the early 1980s (2, 3), the mechanisms by which specific HPV oncoproteins disrupt cellular pathways have been extensively elucidated. HPV types are classified as low-risk or high-risk based on their oncogenic potential (4). While HPV-16 and HPV-18 have been extensively studied, recent epidemiological and molecular studies have underscored the significance of HPV-31 and HPV-52 in cervical cancer incidence, particularly in East Asia and specific regions of Europe (5–8). However, the structural and functional characteristics of the E6 and E7 proteins of HPV-31 and HPV-52 remain poorly characterized.

The oncogenic potential of HPV largely depends on its early proteins, E6 and E7, which facilitate malignant transformation by targeting tumor suppressor pathways (9, 10). E6 binds the p53 tumor suppressor, promoting ubiquitin-mediated degradation and inhibiting apoptosis, while E7 disrupts the retinoblastoma (Rb) pathway to release E2F transcription factors and deregulate cell cycle progression (11–14). Although these mechanisms are conserved among high-risk HPV types, sequence variations in E6 and E7 can lead to differential binding affinities and functional outcomes (15). Recent structural studies have begun to resolve the atomic-level details of HPV-31 and HPV-52 E6 and E7, revealing subtype-specific conformational features that may influence oncogenic potency (16, 18, 19). Nevertheless, a gap remains in the comprehensive bioinformatics characterization of the E6 and E7 proteins of HPV-31 and HPV-52, particularly regarding antigenic epitope prediction—an essential step in vaccine design.

Advances in high-throughput sequencing and computational biology have enabled multidimensional bioinformatics analyses of HPV oncoproteins (16–19). Specifically, homology modeling, molecular docking, epitope mapping, and phylogenetic profiling have uncovered key insights into structural motifs and functional domains of E6 and E7. For instance, Conrady et al. resolved the HPV-31 E6 crystal structure and characterized its interactions with E6AP and p53 (19), whereas Ferenczi et al. conducted phylogenetic and functional analyses of HPV-31 E6 and E7 variants (18). Recent work by Kogure et al. revealed significant intra-patient genomic variability of HPV-31 in cervical cancer and precancer, underscoring the importance of considering viral quasispecies diversity when predicting E6 and E7 epitope profiles (20). Song et al. characterized the genetic variability and phylogeny of HPV-52 E6 and E7 in Sichuan, China, underscoring subtype-specific functional differences relevant to epitope selection (17). Pinheiro et al. conducted a large-scale phylogenomic analysis of HPV-31 across 2,093 genomes, linking specific viral clades to cervical carcinogenesis risk and thereby supporting targeted epitope selection based on subtype phylogeny (21). In summary, prior research has addressed HPV-31 and HPV-52 from various perspectives—sequence diversity (17, 18, 21), structural elucidation (19), and L1 protein-based VLP design (22, 23)—yet none has integrated physicochemical profiling, secondary and tertiary structure modeling, post-translational modification predictions, and B- and T-cell epitope mapping into a single, multilayered framework. Bioinformatics profiling of both subtypes remains incomplete, particularly concerning immunogenic epitope prediction, which is critical for next-generation vaccine design (24).

In this study, the E6 and E7 proteins of HPV-31 and HPV-52 were systematically analyzed using a combination of bioinformatics tools to predict physicochemical properties, post-translational modification sites, secondary and tertiary structures, and to identify potential T-cell and B-cell epitopes. The following hypotheses were tested:

HPV-31 and HPV-52 E6 and E7 proteins exhibit subtype-specific sequence and structural variations that lead to distinct distributions of immunogenic epitopes.The simultaneous application of multiple bioinformatics tools to identical sequences was hypothesized to enhance the accuracy of predicting dominant T-cell and B-cell epitopes in HPV-31 and HPV-52 E6 and E7 proteins.By comparing predicted post-translational modification (PTM) sites with conserved regions, immunogenic regions that may be cross-reactive between subtypes were expected to be uncovered.

Further, it was hypothesized that structural disparities between HPV-31 and HPV-52 E6 and E7 proteins correlate with unique antigenic epitope landscapes, thereby informing the design of future peptide-based vaccines.

Materials and methods

2

Amino acid sequence

2.1

The complete sequence of E6 and E7 oncoproteins of HPV-31 and HPV-52 was available from the National Center for Biotechnology Information(NCBI) database (accession numbers: HPV31 E6 [WAB53637], HPV31 E7 [WAB53638], HPV52 E6 [WAB54303], HPV52 E7 [WAB54304]).

Prediction of protein physicochemical parameters

2.2

Rationale for tool selection and distinctions

2.2.1

To assess basic physicochemical properties of HPV-31/52 E6/E7 proteins, we employed two ExPASy tools:

ProtParam (ExPASy ProtParam v2023.1): We used ProtParam to compute molecular weight, theoretical isoelectric point (pI), extinction coefficient, instability index (II), aliphatic index, and GRAVY (grand average of hydropathicity) in a single run. ProtParam is widely used in viral protein studies because its predictions correlate well with experimentally determined parameters. The instability index (II) quantifies the likelihood of a protein’s stability in vitro, where a value of II > 40 indicates predicted instability (25).

ProtScale (ExPASy ProtScale v2023.1): While ProtParam provides global physicochemical metrics, ProtScale generates residue-level hydrophobicity (Kyte–Doolittle) and hydrophilicity (Hopp–Woods) plots, allowing us to identify local peaks or valleys that may correspond to linear B-cell epitopes. ProtScale employs a sliding-window approach (window size = 7) to generate a continuous hydropathy profile, which ProtParam does not offer (26).

Procedure and statistical processing

2.2.2

The ProtParam calculations were performed in triplicate, and the reported values represent the mean ± standard deviation (SD) of three independent runs.

For the ProtScale analysis, the window width was set to 7 with a default threshold of 0.5. We identified the top three hydrophilicity peaks (using the Hopp–Woods scale) and the deepest hydrophobic valleys (using the Kyte–Doolittle scale) for each protein.

No statistical tests, such as t-tests or ANOVA, were applied because this study is purely predictive, without experimental group comparisons. The results are presented as raw means ± SD for ProtParam values and qualitative hydropathy profiles for ProtScale.

Post-translational modification site prediction

2.3

Rationale for tool selection

2.3.1

NetPhos 3.1 (threshold 0.5): A neural-network–based tool that predicts Ser/Thr/Tyr phos-phorylation sites. We chose NetPhos because it has been benchmarked on short viral proteins with ≥70% accuracy (27). Compared to other open-source servers (e.g., PhosphoSite), NetPhos offers a user-friendly batch interface and provides clear residue-level confidence scores.

MotifScan v2022 (threshold 0.5): Identifies kinase-specific motifs (CK2, PKC, TK, etc.) by searching against curated motif databases (28). We selected MotifScan because it integrates multiple kinase-motif libraries and is particularly suited for mapping short linear motifs adjacent to known functional domains (e.g., LxxLL, LxCxE).

NetNGlyc 1.0 (threshold 0.5): Predicts N-linked glycosylation sites (N-X-S/T motifs) (29). Although E6/E7 proteins rarely undergo glycosylation, we included NetNGlyc to confirm the absence of glycosylation sites—a negative result that supports the cytosolic/nuclear localization of these oncoproteins.

Procedure and output

2.3.2

NetPhos 3.1

2.3.2.1

Submitted each E6/E7 sequence (single sequence mode), extracted residues with score > 0.5.

MotifScan v2022

2.3.2.2

Used default scoring matrices to detect CK2, PKC, TK motifs; only motifs with score > 0.5 were retained.

NetNGlyc 1.0

2.3.2.3

Confirmed that none of the four proteins contained an N-linked glycosylation motif above threshold 0.5.

Signal peptide and transmembrane helix prediction

2.4

SignalP 4.1 (D-score 0.45): Uses a neural network model to predict signal peptide cleavage sites (30). We chose SignalP 4.1 instead of older versions because it offers improved accuracy for proteins lacking obvious signal partners. Its published D-score threshold of 0.45 is recommended for viral oncoproteins.

TMHMM 2.0 (probability threshold 0.5): Predicts transmembrane helices using a hidden Markov model (30). We used TMHMM to verify that E6/E7 do not contain any transmembrane segments, confirming their expected nuclear/cytoplasmic localization.

Secondary structure prediction

2.5

SOPMA v3.0 predicted secondary structure elements (α-helix, β-sheet, β-turn, and random coil) using the default threshold (8% difference, window width = 17). SOPMA’s reported accuracy for viral proteins is ≥70% (31). Compared to alternatives such as PSIPRED, SOPMA provides a residue‐level map that can be directly aligned with predicted epitope regions.

Tertiary structure prediction

2.6

Phyre2 v2.0 (Protein Homology/analogY Recognition Engine) (32) was used for homology modeling of E6/E7 proteins. It leverages experimentally resolved PDB templates and generates high-confidence models for proteins with known homologues (33). Although AlphaFold v3 (2024) can produce de novo predictions, Phyre2’s reliance on validated templates ensures that our HPV E6/E7 models remain directly comparable to prior structural studies (19, 34). This consistency is crucial for accurately mapping predicted epitopes onto known functional domains.

We accepted templates only if they exhibited ≥ 90% sequence coverage and ≥ 99% confidence. Each E6/E7 sequence was submitted in single-sequence mode. For HPV-31 E6, templates c4gizC (coverage 93%, confidence 100%) were chosen; for HPV-31 E7, template d2ewla1 (coverage 50%, confidence 99.8%) was used; for HPV-52 E6, c4gizC (coverage 94%, confidence 100%); for HPV-52 E7, d2b9da1 (coverage 47%, confidence 99.8%).

Template Selection Rationale:

c4gizC: High sequence identity (≥ 90%) with HPV-31/52 E6 in residues 2–144/2–142, respectively (19, 33).

d2ewla1/d2b9da1: Best available templates for E7 with ≥ 99.8% confidence.

Although AlphaFold v3 could produce end-to-end predictions, Phyre2’s reliance on experimentally validated templates (e.g., c4gizC) provides clear alignment evidence and facilitates comparability with existing HPV structural literature (18, 19, 33, 35).

Sequence homology and phylogenetic analysis

2.7

Clustal X 2.0 was chosen for multiple sequence alignment (MSA) because it provides a graphical user interface and allows manual inspection of alignment gaps and conserved motifs. Although other aligners exist (e.g., MUSCLE), Clustal X is widely cited in HPV research and facilitates identification of conserved blocks (≥70% identity).

MEGA 7.0.20 (Molecular Evolutionary Genetics Analysis) was used to construct a Neighbor-Joining phylogenetic tree with 1,000 bootstrap replicates, providing statistical support for each branch. MEGA’s integrated alignment viewer and tree-editing capabilities streamline the generation of publication‐quality phylograms.

We aligned full-length E6/E7 protein sequences from HPV types 16, 18, 31, 33, 35, 45, 52, 56, 58, and 61 using Clustal X 2.0 (gap open penalty = 10; gap extension = 0.1). Evolutionary trees (Neighbor-Joining method, bootstrap = 1,000) were constructed in MEGA 7.0.20 (v7.0.20) to infer phylogenetic relationships. Conserved regions were identified based on ≥ 70% identity across aligned sequences.

Linear epitope analysis of B cells oncoproteins

2.8

We employed four servers to predict linear B-cell epitopes, then selected overlapping regions as dominant candidates:

ABCpred v2.0 (threshold 0.51; peptide length = 16) uses an artificial neural network trained on known linear epitopes (36). We included ABCpred because it has been validated on viral proteins, achieving ~65.9% accuracy (37).

BepiPred 1.0 (threshold 0.35; window = 20) combines hidden Markov models and propensity scales to predict epitopes with a balanced trade-off between specificity and sensitivity (38).

BCPREDS 1.0 (epitope length = 20; specificity = 75%) uses subsequence kernels to identify linear B-cell epitopes; it excels in reducing false positives among random coil regions (39).

SVMTrip v1.0 (threshold 0.51; peptide length = 20) employs a support vector machine algorithm combined with amino acid pair propensity; it outperforms many single‐algorithm tools in independently benchmarked tests (40).

Each E6/E7 sequence was submitted to all four servers in single‐sequence mode. We recorded all predicted peptide segments that surpassed each server’s threshold. Only peptides predicted by ≥ 2 servers were considered for final selection.

Prediction of T-cell epitopes

2.9

CD4^+^ T cell epitopes were predicted using both SYFPEITHI v1.0 (41) and the IEDB MHC II module (42) with HLA-DRB1*15:01 as the reference allele, selected for its 20% frequency in the Chinese population (43). SYFPEITHI is a motif-based predictor that assigns quantitative scores based on known anchor-residue preferences; peptides scoring ≥ 20 were considered strong binders. The IEDB MHC II module generates consensus predictions by integrating multiple algorithms (e.g., NN-align, SMM-align) and has outperformed standalone tools such as TEPITOPE in benchmark studies; CD4^+^ epitopes with a percentile rank ≤ 10 were deemed strong binders.

CD8^+^ T cell epitopes were predicted using the IEDB MHC I module (NetMHCpan 4.1) with HLA-A11:01 and HLA-A02:01—alleles occurring at 18.0% and 15.3% frequency in Chinese individuals, respectively (43). NetMHCpan 4.1 employs a pan-specific neural network to predict peptide binding across diverse HLA-A and HLA-B alleles, consistently outperforming earlier NetMHC versions, especially for less common alleles; CD8^+^ epitopes with a percentile rank ≤ 1 were classified as strong binders. All alleles were chosen based on high-frequency HLA data in the Chinese population (44, 45). The aforementioned methods and corresponding software are summarized in Table 1.

Results

3

Primary structure of HPV-31 and 52 E6 and E7 proteins

3.1

The complete amino acid sequences retrieved from NCBI (HPV-31 E6: 149 AA; HPV-31 E7: 98 AA; HPV-52 E6: 148 AA; HPV-52 E7: 99 AA) are listed below:

HPV-31 E6 (149 AA):

MFKNPAERPRKLHELSSALEIPYDELRLNCVYCKGQLTETEVLDFAFTDL-TIVYRDDTPYGVCTKCLRFYSKVSEFRWYRYSVYGT TLEKLTNKGICDLLIR-CITCQRPLCPEEKQRHLDKKKRFHNIGGRWTGRCIVCWRRPRTETQV

HPV-31 E7 (98 AA):

MRGETPTLQDYVLDLQPEATDLYCYEQLPDSSDEEDVID-SPAGQAKPDTSNYNIVTFCCQCESTLRLCVQSTQVDIRILQELLMGS F GIVCPNCSTRL

HPV-52 E6 (148 AA):

MFEDPATRPRTLHELCEVLEESVHEIRLQCVQCKKELQRREVYKFLFTDLRIVYR DNNPYGVCIMCLRFLSKISEYRHYQYSLYGKTLEERV RKPLSEITIRCIICQTPLCPEEKERH VNANKRFHNIMGRWTGRCSECWRPRPVTQV

HPV-52 E7 (99 AA):

MRGDKATIKDYILD LQPETTDLHCYEQLGDSSDEEDTD GVDRPDGQAEQATSNYYIVTYCHSCDSTLRLCIHSTATDLRTLQQMLLGTLQVVCPGCAR

The physicochemical parameters of the proteins

3.2

Methods brief

3.2.1

ProtParam v2023.1 was used to compute the length, molecular weight, theoretical pI, instability index (II), aliphatic index, and GRAVY. Each value is the mean ± SD of three independent runs.

ProtScale v2023.1 (window size = 7, threshold = 0.5) was used to generate Hopp–Woods hydrophilicity and Kyte–Doolittle hydrophobicity plots to localize potential B-cell epitopes.

All four proteins have a molecular weight >10 kDa, consistent with the reported immunogenic thresholds (46). Instability indices >40 suggest they are intrinsically unstable, potentially influencing antigen processing (37, 46). Negative GRAVY values classify them as hydrophilic, favoring solubility and surface exposure.

Hydrophilicity/hydrophobicity plots (ProtScale) indicate several predicted hydrophilic peaks in the protein sequences (Figure 1). The physicochemical parameters for all four proteins are summarized in Table 2.

Phosphorylation sites: (A) HPV31 E6 (B) HPV31 E7 (C) HPV52 E6 (D) HPV52 E7.

Post-translational modification and subcellular localization predictions

3.3

Methods brief

3.3.1

NetPhos 3.1 (threshold 0.5) was used to predict Ser/Thr/Tyr phosphorylation sites.

MotifScan v2022 (threshold 0.5) was used to identify CK2, PKC, and tyrosine kinase (TK) motifs.

NetNGlyc 1.0 (threshold 0.5) was used to examine possible N-glycosylation sites.

SignalP 4.1 (D-score 0.45) and TMHMM 2.0 (probability 0.5) were used to check for signal peptides and transmembrane helices.

Key findings

3.3.2

The post-translational modification sites and membrane localization of the four proteins are summarized in Table 3. Both E6 proteins have Ser/Thr phosphorylation sites clustered around LxxLL motifs (e.g., S82), suggesting potential regulation of E6AP/p53 binding.

E7 proteins of both subtypes have CK2 sites near the LxCxE Rb-binding motif, suggesting modulation of Rb interaction.

No N-glycosylation, signal peptides, or transmembrane helices were predicted for any of the four proteins, consistent with their known nuclear/cytosolic localization (Figures 2, 3).

Phosphorylation sites: (A) HPV31 E6 (B) HPV31 E7 (C) HPV52 E6 (D) HPV52 E7.

TMHMM analyzed the transmembrane domain of the proteins. (A) HPV31 E6 (B) HPV31 E7 (C) HPV52 E6 (D) HPV52 E7.

Secondary structure predictions

3.4

Methods brief

3.4.1

SOPMA v3.0 (window size = 17, threshold = 8%) was used to determine the percentages of α-helix, β-sheet, β-turn, and random coil.

According to the spatial characteristics of secondary structure, α-helix and β-sheet are not easily disrupted due to hydrogen bonding and are mostly located in the interior of the protein, making them less suitable as antigen-recognizing sites. In contrast, β-turns and irregular curls are primarily protruding structures on the protein surface (47). The specific details of the secondary structures of the four proteins are presented in Table 4. The secondary structure of the HPV-31 E6 protein was analyzed online using SOPMA (Figure 4A). The analysis showed that α-helix accounted for 49.66%, β-sheet for 14.56%, β-turn for 4.43%, and irregular curl for 35.44%. The results indicated that the HPV-31 E6 protein structure is relatively compact (34).

Secondary structure prediction: (A) HPV31 E6 oncoprotein; (B) HPV31 E7 oncoprotein; (C) HPV52 E6 oncoprotein; (D) HPV52 E7 oncoprotein.

The results for the HPV-31 E7 protein showed that α-helix accounted for 25.51%, β-sheet for 22.45%, β-turn for 0%, and irregular curl for 52.04%, as shown in Figure 4B. The results indicated that the HPV-31 E7 protein structure is relatively loose.

For the HPV-52 E6 protein (Figure 4C), α-helix accounted for 54.05%, β-sheet for 10.81%, β-turn for 1.35%, and irregular curl for 33.78%, indicating that the protein structure is relatively compact.

For the HPV-52 E7 protein (Figure 4D), α-helix accounted for 27.27%, β-sheet for 21.21%, β-turn for 0%, and irregular curl for 51.52%, indicating that the protein structure is relatively loose.

Tertiary structure prediction (Phyre2 v2.0)

3.5

Based on Phyre2 outputs (33), high-confidence homology models were obtained for all four proteins (confidence ≥ 99.8%) (Figures 5A–D).

Tertiary structure prediction. (A) HPV31 E6 protein; (B) HPV31 E7 protein; (C) HPV52 E6 protein; (D) HPV52 E7 protein.

HPV-31 E6: The model is based on c4gizC (93% coverage, 100% confidence) (Figure 5A).

HPV-31 E7: The model is based on d2ewla1 (50% coverage, 99.8% confidence) (Figure 5B).

HPV-52 E6: The model is based on c4gizC (94% coverage, 100% confidence) (Figure 5C).

HPV-52 E7: The model is based on d2b9da1 (47% coverage, 99.8% confidence) (Figure 5D).

Key findings

3.5.1

E6 proteins are helix-rich and compact, with fewer β-turns, suggesting that most linear epitopes lie in random coil loops.

E7 proteins contain ≥ 50% random coil, indicating extensive surface exposure and many potential linear epitopes.

HPV-31 and HPV-52 E6/E7 structures are highly conserved overall, with only minor local deviations that may underlie subtype-specific immunogenic differences.

Homology and phylogenetic analysis (Clustal X 2.0 & MEGA 7.0)

3.6

Amino acid identity and conserved regions

3.6.1

Multiple sequence alignment of E6 proteins (HPV-16, 18, 31, 33, 35, 45, 52, 56, 58, 61) revealed conserved motifs at positions 8–15, 25–34, 41–77, 79–89, 96–112, 114–141 for HPV-31 E6, and 8–16, 25–31, 41–56, 59–69, 71–79, 81–89, 101–107, 109–119, 123–125, 130–136 for HPV-52 E6 (Figure 6A). E7 proteins exhibited conserved regions at 1–17, 20–28, 30–36, 38–45, 52–77, 82–87, 89–94 (HPV-31) and 10–15, 24–28, 30–36, 39–46, 53–59, 62–70, 76–96 (HPV-52) (Figure 6C).

Homology and molecular evolution analysis. (A) Homology analysis of E6 proteins of HPV; (B) The molecular evolutionary tree of E6 proteins of HPV; (C) Homology analysis of E7 proteins of HPV; (D) The molecular evolutionary tree of E7 proteins of HPV.

Conserved regions overlap predicted epitope regions, suggesting potential cross-reactivity among related types (48). The HPV−31 E6 45–53 region aligns with the HPV−16 E6 45–53 region, indicating possible shared immune responses.

Phylogenetic tree construction

3.6.2

Neighbor-Joining trees (bootstrap = 1,000) placed HPV-31 E6 in a close clade with HPV-35 E6 (Figure 6B), and HPV-52 E6 in a close clade with HPV-33 E6. For E7, HPV-31 clustered with HPV-16, while HPV-52 clustered with HPV-33 (Figure 6D).

Linear epitopes of B cells

3.7

Methods brief

3.7.1

Tools: ABCpred v2.0 (peptide length = 16; threshold = 0.51), BepiPred 1.0 (threshold = 0.35), BCPREDS 1.0 (peptide length = 20; specificity = 75%), and SVMTrip v1.0 (peptide length = 20; threshold = 0.51).

Criterion: Retain only peptides predicted by ≥2 algorithms and restrict to loop/turn regions identified by SOPMA.

After excluding α-helix and β-sheet regions, the top five predicted epitopes per method were compared. Using the four B-cell prediction tools, overlapping epitopes (predicted by ≥2 servers) were identified as dominant (Supplementary Tables 1–16). After cross-referencing, the dominant B-cell epitopes were Table 5:

HPV-31 E6: 55–61 (RDDTPYG), 112–116 (PEEKQ), 125–131 (FHNIGGR)HPV-31 E7: 8–17 (LQDYVLDLQPEATDLYC), 16–20 (QPEAT), 29–41 (PDSSDEEDVIDEP), 42–48 (AGQAKPDT)HPV-52 E6: 110–119 (LCPEEKERHV), 129–141 (MGRWTGRCSECWR)HPV-52 E7: 11–19 (YILDLQPET), 23–27 (HCYEQ), 29–38 (GDSSDEEDTD), 36–48 (DTDGVDRPDGQAE)

Key findings

3.7.2

HPV-31 E6 candidate epitopes (e.g., 55–61 RDDTPYG) are located in a random coil adjacent to LxxLL, suggesting potential for neutralizing antibodies.

The HPV-31 E7 region 29–41 (PDSSDEEDVIDEP) is consistently predicted by four methods and is located within a highly exposed coil loop.

The C-terminal loops of HPV-52 E6/E7 (e.g., 129–141 in E6, 36–48 in E7) are strong candidates for B-cell epitopes.

Linear epitopes of T cells

3.8

CD4+ T cell epitope prediction (HLA-DRB1*1501)

3.8.1

The SYFPEITHI and IEDB MHC II tools (percentile rank ≤ 10; positive control) were used. Supplementary Tables 17–20 present the top five predictions. The final dominant CD4^+^ epitopes (overlapping high-scoring predictions) are as follows:

• HPV-31 E6: 45–53 (FAFTDLTIV), 72–80 (KVSEFRWYR).- HPV-31 E7: 7–15 (TLQDYVLDL), 11–19 (YVLDLQPEA), 82–90 (LLMGSFGIV).- HPV-52 E6: 45–53 (FLFTDLRIV), 82–87 (SLYGKT).- HPV-52 E7: 84–90 (MLLGTLQ), 53–59 (NYYIVTY), 11–19 (YILDLQPET).

CD8+ T−cell epitope prediction (HLA-A1101, A0201)

3.8.2

IEDB MHC I binding (NetMHCpan 4.1; percentile rank ≤ 1) was used. Supplementary Tables 21–24 present the results. The final dominant CD8^+^ epitopes are as follows (Table 6):

• HPV-31 E6: 82–90 (SVYGTTLEK; HLA-A1101 rank 0.01), 45–53 (FAFTDLTIV; HLA-A0201 rank 0.93)- HPV-31 E7: 7–15 (TLQDYVLDL; HLA-A0201 rank 0.09), 37–46 (VIDSPAGQAK; HLA-A1101 rank 0.33)- HPV-52 E6: 86–94 (KTLEERVRK; HLA-A1101 rank 0.01), 18–26 (VLEESVHEI; HLA-A0201 rank 0.03)- HPV-52 E7: 84–92 (MLLGTLQVV; HLA-A0201 rank 0.08), 51–59 (TSNYYIVTY; HLA-A1101 rank 0.74)

Notably, the overlapping T-cell epitope 45–53 appears in both E6 proteins and is conserved between HPV-31 and HPV-52, suggesting a promiscuous HLA-binding region that could elicit cross-type T-cell responses.

Discussion

4

In this study, integrative bioinformatics approaches were employed to analyze the E6 and E7 proteins of HPV-31 and HPV-52, identifying key structural features and dominant antigenic epitopes. The key findings and their biological implications are addressed in the subsequent sections.

Physicochemical properties and implications for immunogenicity

4.1

Viral proteins with molecular weights exceeding 10 kDa typically exhibit sufficient immunogenicity for epitope recognition (46, 49). All four E6 and E7 proteins of HPV-31 and HPV-52 exceed this threshold (17.8–18.0 kDa) and are classified by ProtParam as “unstable” (instability index > 40), a feature associated with increased post-translational susceptibility and potential antigenicity (37, 50, 51). Negative GRAVY scores categorize these proteins as hydrophilic, thereby promoting solubility and enhancing epitope exposure (52). These properties correlate with an enhanced potential for antigen presentation, which is critical for vaccine design.

Post-translational modifications and functional context

4.2

Predicted phosphorylation sites were mapped to residues involved in the interactions of E6 and E7 with host regulators. For instance, conserved serine residues (S82 in both E6 proteins) reside within the LxxLL-binding pocket, which is crucial for E6AP-mediated p53 degradation (15, 19). CK2 phosphorylation motifs overlapping this region may modulate binding affinity and subsequent ubiquitination (10, 17). Similarly, E7 CK2 sites (e.g., residues 7–10 encompassing the LxCxE motif) likely regulate Rb binding, contributing to cell cycle dysregulation (11, 13). PKC sites adjacent to the C-terminal zinc-finger (E6 133–135) may influence nuclear localization and stability (15). These in silico insights align with experimental evidence showing that kinase-mediated phosphorylation directly alters oncoprotein function (10, 19).

Secondary/tertiary structures and template selection

4.3

SOPMA analysis reveals that the E6 proteins are predominantly composed of α-helices (49.66% in HPV-31; 54.05% in HPV-52), suggesting compact cores that may shield specific epitopes. In contrast, the E7 proteins exhibit a higher proportion of random coils (52.04% and 51.52%, respectively), indicating flexible surface regions conducive to antibody binding (37, 53). Previous studies have shown that random coils frequently coincide with B-cell epitope hotspots (52, 53), supporting our predictions of dominant linear B-cell epitopes within coil-rich segments, such as residues 8–17 (HPV-31 E7) and 23–27 (HPV-52 E7).

Homology models generated by Phyre2 (confidence > 99.8%) confirm conserved structural motifs, including zinc-binding Cys motifs, consistent with experimental structures (19, 40). The 3D models generated by Phyre2, validated by high confidence scores, display conserved zinc-finger motifs and binding pockets. While AlphaFold3 (2025 release) could generate full-length models, Phyre2’s template-based approach allowed for a direct comparison with known E6/E7 structures. We selected Phyre2 templates (c4gizC/d2ewla1/d2b9da1) due to their high sequence identity (>50%) and prior experimental validation (18, 19).

Homology and evolutionary insights

4.4

Multiple sequence alignment and phylogenetic analysis position HPV-31 E6 closely with HPV-35, and HPV-52 E6 with HPV-33, while E7 clusters similarly with HPV-16 and HPV-33 (18, 24). Conserved regions (e.g., E6 positions 41–77; E7 positions 52–77) overlap with predicted T-cell epitopes, suggesting potential cross-reactivity and cross-protection among high-risk HPV types (16, 17). This cross-immunity is essential for the design of multivalent vaccines targeting broad high-risk HPV coverage (6).

Antigenic epitope identification and validation potential

4.5

Dominant B-cell epitopes were identified (e.g., HPV-31 E6: 55–61, 112–116, 125–131; HPV-52 E7: 23–27, 29–38, 36–48) and T-cell epitopes (e.g., HPV-31 E6: 45–53; HPV-52 E6: 86–94), predicted by multiple algorithms (ABCpred, BepiPred 1.0, BCPREDS, SVMTriP) (37, 54, 55). CD8^+^ epitopes, such as HPV-31 E7: 7–15 (TLQDYVLDL), exhibited a strong binding affinity to HLA-A0201 (IEDB rank 0.09), consistent with known CTL responses against HPV-16 E7 (11, 56). CD4^+^ epitopes (e.g., HPV-52 E7: 11–19) exhibited favorable binding to HLA-DRB1*1501, which is crucial for helper T-cell activation (42). These in silico predictions align with experimental data linking epitope immunodominance to surface accessibility and structural features (54, 55). Subsequent empirical validation, such as peptide-MHC binding assays and T-cell activation studies, is necessary (42, 56).

Comparison with previous studies

4.6

Previous studies have characterized the sequence variability of HPV-31/52 (17, 18, 21) and resolved individual E6 crystal structures (19). Kogure et al. further demonstrated that HPV-31 genomes exhibit significant intra-patient heterogeneity (20), suggesting that E6 and E7 epitopes may evolve during disease progression. However, to date, no study has integrated physicochemical properties, post-translational modification site prediction, secondary and tertiary structure modeling, and multilayered immunoinformatic epitope mapping for both E6 and E7 of HPV-31 and HPV-52 into a single comprehensive analysis. Our work addresses this gap by correlating predicted phosphorylation sites with functional motifs (e.g., LxxLL, LxCxE) (27, 57) and mapping B- and T-cell epitopes to conserved, surface-exposed regions identified through structural modeling. Furthermore, Song et al. and Firdaus et al. have highlighted the immunogenic potential of HPV-52 (17, 22, 23), particularly in Asian populations, thus validating the public health relevance of our subtype-specific epitope predictions. Kesheh et al. proposed region-tailored multivalent vaccine designs based on L1 gene diversity (58), offering translational context for our E6 and E7-based epitope candidates.

Application to vaccine design

4.7

Although this study did not experimentally construct virus-like particles (VLPs) or multivalent peptide vaccines, the predicted epitopes provide a foundation for rational vaccine design:

Cross-subtype conserved CD8+ epitopes

4.7.1

The E6 45–53 segment in HPV-31 (FAFTDLTIV) and HPV-52 (FLFTDLRIV) exhibits strong binding affinity for HLA-A0201 and HLA-A1101 (IEDB percentile ≤ 1) and is highly conserved across high-risk types, making it an ideal candidate for inclusion as a universal cytotoxic T-lymphocyte (CTL) epitope in a multi-epitope DNA or peptide vaccine.

Helper T-cell (CD4+) epitopes

4.7.2

E6 72–80 (KVSEFRWYR) in HPV-31 and E6 82–87 (SLYGKT) in HPV-52 exhibit moderate binding affinity to HLA-DRB1*1501 (IEDB percentile ≤ 10) and could be fused with CTL epitopes into a single recombinant protein or synthetic long peptide construct to enhance helper T-cell responses, as suggested by He et al (57).

B-cell neutralizing epitopes on VLP platforms

4.7.3

The B-cell epitope HPV-31 E6 55–61 (RDDTPYG) and HPV-52 E6 110–119 (LCPEEKERHV) reside in exposed random coil regions. Firdaus et al. successfully inserted analogous linear epitopes into the L1 VLP platform to elicit neutralizing antibodies (22), supporting the strategy of grafting these peptides onto L1 VLPs to generate subtype-specific antibody responses.

Multivalent peptide/protein vaccine constructs

4.7.4

Building on Firdaus et al.’s reverse vaccinology design for HPV-52 L1 (23), one could concatenate top CD4^+^ and CD8^+^ epitopes (e.g., E6 45–53, 72–80; E7 7–15) with appropriate linkers and trafficking signals to create a chimeric protein capable of eliciting robust humoral and cell-mediated immunity in preclinical HLA-transgenic mouse models.

Limitations and future directions

5

Although the integrated in silico pipeline provides a comprehensive epitope landscape, experimental validation—such as peptide-MHC binding assays, ELISpot, and crystallographic studies—is crucial to confirm immunogenicity (54, 55). Additionally, molecular dynamics simulations could refine epitope conformations and assess stability within MHC binding grooves (32, 51). This study relies solely on in silico predictions and lacks direct in vitro or in vivo validation, representing a primary limitation. Pinheiro et al. confirmed that certain E6 and E7 regions correlate with cervical cancer aggressiveness at the genomic level (21), yet these findings require empirical confirmation through immunological assays. Kogure et al. observed intra-patient HPV-31 variants across different lesion stages (20), emphasizing the need to validate epitope immunogenicity across clinical time points. Future studies should involve:

Experimental binding assays

5.1

Use ELISPOT or flow cytometry with peptide-stimulated peripheral blood mononuclear cells (PBMCs) from HLA-typed donors to validate CD4^+^ and CD8^+^ T-cell responses against the predicted epitopes.

Antibody neutralization studies

5.2

Synthesize candidate B-cell epitopes (e.g., HPV-31 E6 55–61; HPV-52 E6 110–119) and assess their ability to induce neutralizing antibodies in ELISA or pseudovirus neutralization assays.

Animal model validation

5.3

Evaluate peptide-based or VLP-based vaccine constructs (e.g., insertion of linear epitopes into L1 VLPs, as demonstrated by Firdaus et al., 2023) in HLA-transgenic mouse models to measure protective efficacy against HPV-induced tumorigenesis.

In summary, the integrative bioinformatics analysis illuminates subtype-specific structural and immunogenic features of HPV-31 and HPV-52 E6 and E7 proteins, laying the groundwork for experimental validation and rational vaccine design aimed at reducing the HPV-associated cervical cancer burden.