EpiCurator: an immunoinformatic workflow to predict and prioritize SARS-CoV-2 epitopes

Genome retrieval and protein annotation

Genome sequences of 1,652 SARS-CoV-2 genome isolated in Brazil were retrieved from the GISAID database (Elbe & Buckland-Merrett, 2017) available from December 2020 to April 2021 (Table S1). The lineage distribution of the retrieved genomes includes P.1 (Gamma, n = 770), P.2 (Zeta, n = 525), B.1.1.28 (n = 223) and B.1.1.33 (n = 136). These genomes were analyzed by The Viral Annotation Pipeline and iDentification (VAPiD) (Shean et al., 2019) to determine the amino acid sequence for all SARS-CoV-2 proteins using as reference NC_045512.2 from the NCBI database (https://www.ncbi.nlm.nih.gov/). The final genome processing includes clustering the amino acid sequences for each protein using the CD-HIT package (Li & Godzik, 2006) with a 100% sequence identity threshold.

Spike protein comparison

Clustal Omega and MView tool performed a homology analysis for the Spike protein (Sievers & Higgins, 2018); (Brown, Leroy & Sander, 1998) among 100 random samples of Brazilian lineages and the VOCs, retrieved from the GISAID database (Elbe & Buckland-Merrett, 2017), to measure the identity of the sequences to Gamma (P.1) samples.

Epitope prediction

The T and B cell epitopes prediction was performed using the softwares NetCTL, NetMHCpan, and NetMHCIIpan (https://services.healthtech.dtu.dk/) that allow the high-throughput computing analysis. These predictors support FASTA files containing amino acids marked with “X” in the sequence expanding the possibility of public genomes analysis, even with minor sequencing errors. This prediction features decreasing the pre-processing sequence steps and avoids the false-positive epitopes provided by joining sequences.

Prediction of SARS-CoV-2 epitopes

To predict CD8+ T cell epitopes, the FASTA sequences of SARS-CoV-2 proteins are processed using NetCTL v1.2 (Larsen et al., 2007). First, the sequences and the HLA class I supertypes, provided by the software, are submitted to select 9-mer peptides (Chen et al., 1994; Sakaguchi et al., 1997; Gfeller et al., 2018), then peptides with amino acids marked with “X” are removed (Fig. 1A). The predicted peptides in NetCTL are further processed using NetMHCpan v.4.0 software (Hoof et al., 2009; Jurtz et al., 2017) to identify epitopes with strong binding affinity to HLA class I alleles. The prediction parameter is based on the predicted percentile rank ≤ 0.5% and half-maximal inhibitory concentration (IC50) < 500 nM (Chen et al., 1994; Sakaguchi et al., 1997; Gfeller et al., 2018) (Fig. 1A). To assess binding affinity, alleles of HLA-A, B, and C loci were selected from Allele Frequency Net Database (Gonzalez-Galarza et al., 2020) by their Brazilian population frequency > 5% (Table S2).

To predict CD4+ T cell epitopes and estimate binding affinity to HLA class II molecules, the FASTA sequences are processed using NetMHCIIpan v3.2 software (Greenbaum et al., 2011). This tool selects the epitopes with 15-mer peptide lengths based on the predicted percentile rank ≤ 2.0% and IC50 < 500 nM (Fig. 1A). For the affinity prediction, the loci HLA-DRB1, HLA-DPA1-DPB1, and HLA-DQA1-DQB1 were selected by the phenotypic frequency > 5% in the Brazilian population (Table S2) from Allele Frequency Net Database (Gonzalez-Galarza et al., 2020). For both cells (CD8+ and CD4+ T cell), the epitopes with broad HLA affinity coverage (≥3 alleles) are filtered (Fig. 1A).

The prediction of linear B cell epitopes from SARS-CoV-2 structural proteins was performed by BepiPred v.2.0 software (Jespersen et al., 2017), with a threshold of 0.5 (corresponding specificity > 0.817 and sensitivity < 0.292) (Fig. 1A). Only the epitopes with more than seven and less than 50 amino acid residues in length and sequence without amino acids marked with “X” are considered for subsequent curation analysis.

Accurate selection of epitopes - EpiCurator

We developed a rigorous analysis workflow for accurate selection of epitopes, the EpiCurator, bringing together a set of filters according to pre-established criteria, optimizing the analysis since it brings tools available to use by high-throughput computing architectures. They also group a series of analyses to guarantee the selection of unpublished and qualified epitopes. The workflow allows the identification of epitopes by the following analysis: conservancy, homology with the human genome, the overlap between epitopes of different classes, and the identification of epitopes previously published in PubMed Central^® (PMC) or available in the IEDB database (Vita et al., 2019) beyond identifying their protein coordinates and mutations (Fig. 1B).

Epitope properties

Genomes for EpiCurator pairwise comparison validation

Protein sequences for SARS-CoV-2 isolates reported by Crooke et al. (2020) were identified and retrieved from the Virus Pathogen Resource (ViPR) database (n = 641, 635); additionally, six genome sequences reported by Kiyotani et al. (2020), two sequences for N and S protein reported by Chen et al. (2020), and five S sequences reported by Chukwudozie et al. (2021) were retrieved from NCBI GenBank. These samples were processed using prediction parameters of the HLA alleles reported by the authors with the approach reported in this Methods section (Epitope prediction and Accurate selection analysis). Regarding item Accurate selection of epitopes - EpiCurator of the Methods, we used only the three main analysis steps (conservancy, human homology and IEDB matching) to compare the selected epitopes ensemble by the papers, and ones chosen independently by our approach, leaving out the EpiMiner analysis since all the epitopes are published.

Multi-epitope construct and structural modelling for EpiCurator validation

The RBD epitopes (n = 11) were used to construct a multi-epitope sequence connected by specific linkers (Fig. S1). The linker aimed to separate the epitopes, so that it improved their expression, folding and stability beyond to prevent their fusion and facilitate the immune processing of antigen (Arai et al., 2001; Kar et al., 2020). In the multi-epitope arrangement are also added an adjuvant (Escherichia coli 50S ribosomal protein L7/L12 (UniProt P0A7K2)) in the N-terminal sequence and a histidine hexamer in the C-terminal portion (Fig. 1D, Fig. S1).

Linear epitopes prediction from SARS-CoV-2 genome

Three main epitopes prediction comprise HLA Class I (CD8+ T cell), HLA Class II (CD4+ T cell), and B cell. To predict potential CD8+ T cell epitopes, the NetCTL and NetMHCpan, predictive algorithms were performed for all proteins annotated from 1,652 SARS-CoV-2 genomes (see Methods). Thus 9-mer epitopes were predicted with strong binding affinity assigned by percentile scores (rank) ≤0.5% across HLA class I alleles frequent in the Brazilian population (Fig. 1A). Cumulatively reaching the prediction of 5,261 HLA class I epitopes with HLA promiscuity, binding to ≥3 HLA alleles, implying broad population coverage. The affinity to HLA-C alleles matches the highest epitope binding locus suggesting it could be the best grooves for these epitopes and would lead to the activation of the cell-mediated response (Table 1).

We also sought to predict potential 15-mer epitopes with binding affinity to HLA class II, using the NetMHCIIpan software (Fig. 1A), reaching 7,649 candidate HLA class II epitopes from the SARS-CoV-2 genomes. They have a strong binding affinity to ≥3 HLA alleles across a reference panel of HLA molecules (see Methods). These epitopes showed a preferential affinity for alleles from locus HLA-DRB1 (n = 3, 702), suggesting that HLA-DRB1 alleles could lead these epitopes to activate the cellular immune response (Table 1).

To pairwise the cellular and humoral immune responses activation, the additional prediction of linear B cell epitopes was performed by BepiPred algorithm identifying 257 epitopes from SARS-CoV-2 structural proteins (Table 1).

Quality assessment analysis of the EpiCurator

A robust and refined accurate selection of epitopes is crucial to improve the development of peptide-derived vaccines. To this goal, the EpiCurator brings together six analyses (Conservancy, Human homology, Epitopes overlap, EpiMiner, IEDB matching, and Mutation screening) to the accurately selected epitopes (Fig. 1B).

To the workflow’s quality assessment analysis, we sought to characterize the number of epitopes taken for the main analysis. The epitopes conservancy selected 29.23% of the predicted epitopes with 100% identity across at least 90% of the SARS-CoV-2 samples. Nevertheless, 70.77% of the predicted epitopes do not correspond to the conservancy parameter (Fig. 2, Fig. S2). Human homology analysis identified only 0.1% of the predicted epitopes sequence in the human genome, keeping 99% accurately selected, unmatched with human sequences (Fig. 2, Fig. S2). In addition, the screening by the previously published epitopes (EpiMiner and IEDB matching analysis) allowed selected > 80% of new epitopes, since 17.57% of them have already been described in the literature and/or IEDB server (Fig. 2, Fig. S2). To confirm the effectiveness of the EpiMiner we provide all the articles IDs and DOI in which the epitopes were found in Table S4.

To further assess the quality of the EpiCurator analysis, a pairwise comparison validation was performed. The comparison used NGS data reported by four papers with substantial similarity with our methodologies (Kiyotani et al., 2020; Chen et al., 2020; Crooke et al., 2020; Chukwudozie et al., 2021). The prediction step allows identifying the same HLA class I and HLA class II-restricted T cell epitopes. At the same time, the accurate selection provided by the EpiCurator gathered 51.1% of the paper’s epitopes (Table S5). The conservancy across SARS-CoV-2 genomes samples reached 95,6% of the paper’s epitopes, with only 4.4% not conserved by our parameters. The human homology step certifies the absence of the article’s selected epitopes in the Human genome since the analysis does not recognize the homology of any of them. Interestingly, the IEDB matching step identified 45.6% of the paper’s epitopes as having already been published on the IEDB server showing the importance of this step analysis if we want to describe the epitopes by the first time (Table S5).

Properties of accurately selected epitopes

Assembling the workflow analysis results, 199 (3.78%) HLA class I-restricted T cell epitopes were selected from SARS-CoV-2 proteins, mainly identified in ORF1ab (n = 154 (77.4%)) and Spike glycoprotein (n = 15 (7.5%)) (Table 2, Fig. 3A). The epitopes keep a high affinity for HLA-C alleles (Fig. 3C). The HLA class II epitopes’ accurate selection reached 153 (2%) epitopes, also mainly identified in ORF1ab (n = 111 (72.5%)) and Spike glycoprotein (n = 22 (14.4%)). The HLA class II epitopes showed an affinity prevalence for the HLA-DPA1-DPB1 haplotypes (Fig. 3D). The comprehensive workflow analysis for B cell epitopes selected 14.6% of predicted epitopes with 60% of the epitopes from the Spike glycoprotein (Table 2). Interestingly, all epitopes have high antigenicity (0.92 ± 0.30), are non-toxic and non-allergenic (Table S3).

Furthermore, the T cell epitopes immunogenic potential was assessed to characterize the capability of inducing in silico protective immune responses. The analysis identified the profile of HLA class I and HLA class II epitopes respectively as follow: IL-4 inducer activity (26.6% and 23.8%), IFNγ production (22.1% and 39.7%), immunomodulatory activity (2% and 3.3%), and proinflammatory activity (78.9% and 81.4%) (Table S3).

In addition, to assess the capability to be an in silico Brazilian epitopes candidate, we sought to estimate their distribution among the Brazilian circulating lineages at the time of analysis (P.1, P.2 and B.1.1.28). Firstly, we identify the number of epitopes in common across P.1, P.2 and B.1.1.28 lineages identifying the highest proportion of HLA class II epitopes (61.4%), following by HLA class I (46.7%) and B cell (14.3%) (Fig. 3B). Thereafter, considering the most representative lineage in Brazil, at the time data was retrieved (P.1), remarkably, proportions highest than 50% for all the epitopes were identified, with 70.1% of HLA class II epitopes, 69.8% of HLA class I, and 51.4% of B cell epitopes in P.1 (Fig. 3B).

In the last analysis, considering the HLA promiscuity epitopes selection, the population coverage for the T cell epitopes associated with their respective HLA allele binding (Table S3) was estimated by the IEDB server. Notably, the T cell epitopes have a wide population coverage, presenting 99.52% (HLA class I epitopes) and 100% (HLA class II epitopes) of cumulative Brazilian population coverage and 98.45% (HLA class I epitopes) and 99.87% (HLA class II epitopes) of worldwide population coverage (Table S6).

Epitope-specific RBD Spike as a baseline for validation of EpiCurator selection

The EpiCurator allows the selection of the majority of the epitopes in the ORF1ab and Spike glycoprotein. The SARS-CoV-2 Spike glycoprotein has greater prominence concerning the virus and host interaction and has been the main target in epitopes prediction (Shang et al., 2020; Walls et al., 2020). Therefore, we assumed the spike epitopes as a baseline to validate the EpiCurator accurate selection. The Spike epitopes (n = 58) were mainly in the N-terminal domain (NTD) (53.4%) and receptor-binding domain (RBD) (18.9%) (Fig. 4).

The RBD epitopes (n = 11) comprise three of HLA class I (461-NYNYRYRLF-469, 506-QSYGFQPTY-514, 516-FGYQPYRVV-524), three of HLA class II (322-EKGIYQTSNFRVQPI-336, 322-EKGIYQTSNFRVQPR-336, 443-TGCVIAWNSKNLDSK-457), and five B cell epitopes (332-RVQPTES-338, 424-APGQTGK-430, 424-APGQTGT-430, 455-DSKVGGNYN-463, 474-LKPFERD-480) (Fig. 4).

To evaluate the robust and refined RBD epitopes selection, we perform a diversified analysis. They unveil the highest antigenicity (1.1 ± 0.27) (Fig. 5A) reflecting their ability to bind molecules of adaptive immunity. Notably, their affinity for several HLA alleles are high (0.49 ± 0.15) (Fig. 5A) and achieves promiscuity with ≥ 5 HLA allele binding per epitope (Fig. 5B). The RBD epitopes HLA data pairwise with population coverage that remains at around 80% in the Brazilian and worldwide population (Fig. 5B). To further characterize the RBD epitopes, they reach the conservancy of over 99% across SARS-CoV-2 genomes samples of lineages P.1, P.2 and B.1.1.28 published on GISAID (See Methods) (Table S3).

The high conservancy observed suggests that epitope sequences are shared in the circulating Brazilian lineages. Considering the natural homology among the SARS-CoV-2 lineages but regardless of the epitopes’ conservancy with other lineages, we expand the identification of the RBD epitopes for diversified samples available in GISAID at the time of writing. Interestingly these epitopes were identified in more than 1 million samples for around 1 thousand lineages of GISAID (Table S7). These findings include some samples of current VOCs: Alpha (B.1.1.7), Beta (B.1.351), and Delta (B.1.617.2) (Table S7). Concern about the natural lineages’ homology mentioned, we perform a comparison among spike glycoprotein of described Brazilian lineages and VOCs. The evaluation of random 100 samples available in GISAID of each lineage shows 99.1% of identity with 99.9% of coverage confirming the sequence lineages similarity (Table S8).

One of the mentioned properties of RBD epitopes is the high affinity for HLA alleles. To validate these findings the in silico structural binding performance was assessed by the PepDock docking tool. Thus, the RBD epitopes and the most major HLA alleles in studies with COVID-19 (see Methods) were structurally bonded in complex HLA-epitope (Fig. S3). All complexes were seen with high-performance scores demonstrating in silico epitope effectiveness in activating the cell-mediated immune response via MHC presentation. Their high estimated accuracy selected the two most significant HLA-epitope complexes (Table S9, Fig. 6). The HLA class I epitope 516-FGYQPYRVV-524 in complex with HLA-A*01:01, HLA-A*02:01, HLA-B*08:01, and HLA-C*12:03; and the HLA class II epitope 443-TGCVIAWNSKNLDSK-457 in complex with HLA-DRB1*03:01 and the HLA-DRB1*12:02 alleles (Table S9, Fig. 6). Additionally, the best measured free energy of the complexes was ΔG = −23 kcal/mol for 461-NYNYRYRLF-469 HLA class I epitope in complex with HLA-A*01:01 and ΔG = −31.6 kcal/mol for 322-EKGIYQTSNFRVQPR-336 HLA-class II epitope in complex with the HLA-DRB1*04:01 allele (Table S9, Fig. 6). An additional analysis was performed to assess the specific residues involved in the structural binding allowing to identify the specific amino acids in contact between the epitopes and the HLA alleles (Table S9).

Multi-epitope construct for in silico validation of RBD epitopes

To further validate the RBD epitopes and consequently confirm the effectiveness of EpiCurator in providing a curated selection, the epitopes were inserted in a multi-epitope construct. It consists of 287 amino residues, including 11 RBD selected epitopes joined by linkers (Fig. S1). The structural appraisal of the secondary structure predicted by the PSIPRED server revealed 31.7% alpha-helix, 13.5% beta-strand and the disordered region was 12% (Fig. S4). Some evaluation of the linear sequence of multi-epitope were identified as follows: high antigenicity by the Vaxijen 2.0 (0.87) and the AntigenPro (0.93), good solubility by SolPro (0.91), Protein-sol (0.66), and SOLart (65.80%), non-allergenic feature by AllerTOP and AllergenFP, hydrophilic property (hydropathicity = −0.438) and stability (instability index = 20.56), with pI 7.03, and molecular weight calculated to be 29 kDa.

Another assessment of linear sequence includes the in silico immunogenic profile provided by the C-IMMSIM immune server. The analysis showed a gradual increase of IgM, IFN-γ (associated with both CD8+ T cell and CD4+ Th1 response as shown in Fig. 7A), and IL-2 level after each multi-epitope exposure indicated an elevated immune response (Fig. 7B). Besides, an adequate generation of both IgG1 and IgG2 was shown (Fig. 7C) with high levels of clonal proliferation of B cell and T cell population (Figs. S5A and S5B). In addition, the development of in silico immune memory was assessed by the abundance of different types of B cells and T cells (Figs. S5C and S5D).

The favorable in silico immunogenic profile of linear construct led to a multi-epitope 3D modelling to assess TLR binding capability adding validation for the curated selection provided by the EpiCurator. In this context, we modeled and evaluated the 3D multi-epitope to choose the best model for the TLR binding assay. The best predicted tertiary structure returned by RaptorX has an RMSD score of 11.218 (Fig. S6A), 62% of residues are predicted to be exposed and 30% are predicted to be disordered (Fig. S6B). The refinement by GalaxyRefine allows the choice of the best model based on its quality scores available in Table S10. This model highlights the epitopes according to their type (B cell, HLA class I T cell and HLA class II T cell) (Fig. S7). The overall model quality of the refined structure was further validated in the ProSA web (Fig. S8A), showing the energy results by position (Fig. S8B) with good local quality since all energy values of the residues are negative.

Considering the accuracy of the refined model, a docking between our multi-epitope and the immune receptors (TLR3 and TLR4) was performed to check stability and binding affinity. The best complexes (Fig. 8) for each receptor had free energy values of ΔG = −1,048 kcal/mol (TLR3) and ΔG = −1020 kcal/mol (TLR4), indicating a high binding affinity. In addition, the stability and physical movements of the complexes were confirmed using molecular dynamics simulation in the iMODS server. The complex results are presented in Figs. S9 and S10.