Epitopes, B-Lymphocyte

Let i = (i₁, i₂, ..., i_w) denote a sequence of amino acids, which has been extracted from a protein sequence. Let j denote the position in this window, j = 1...w. On basis of i, the hidden Markov model predicts if the center position of the window is annotated as part of an epitope. In the N- and C-termini, parts of the extracted windows are exceeding the terminals. For these residues, the character 'X' is used, which does not count when the hidden Markov model is used for the predictions. The prediction score for a window is given by
which is the log odds of the residue at the center position of the window is being part of an epitope (Epitope model) as opposed to if it is occurring by chance (Random model).
To construct the Random model, background frequencies of the Swiss-Prot database [23 (link)], q_i, is used. For the Epitope model, p_i,j is the effective amino acid probability of having amino acid i at position j according to the model.
To calculate the values of p_i,j, all windows, for which their center position is annotated as part of an epitope, are extracted from atraining data set. Again, if an extracted window exceeds the N or C terminal, the character 'X' is used, which does not count when calculating the parameters.
These extracted peptide windows form a matrix of aligned peptides of the width w. From this alignment, p_i,j is calculated as the pseudo count corrected probability of occurrence of amino acid i in column j, estimated as in [24 (link)]. To make the pseudo count correction, pseudo count frequencies, g_i,j, are calculated. They are given by
where p_k,j is the observed frequency of amino acid k in column j of the alignment [25 (link)]. The variable b_i,k is the Blosum 62 substitution matrix frequency, e.g. the frequency of which i is aligned to k [26 (link)].
To give an example of using (2), let the window size, w = 1. The model is then only covering residues, which are annotated as being part of linear B-cell epitopes. If the observed peptides consists of the following single amino acid sequences L and V, with the frequencies p_L,1 = 0.5 and p_V,1 = 0.5, then the pseudo-count frequency for e.g. I is given by
The effective amino acid frequencies are calculated as a weighted average of the observed frequency and the pseudo count frequency,
Here, α is the effective number of sequences in the alignment - 1, and β is the pseudo count correction [25 (link)], which is also called the weight on low counts. To finish the calculation example, let β be very large as it is in this work. Then p_I,1 ≈ g_I,1 = 0.14.
Note that we shall use the term hidden Markov model throughout this work to refer to the weight matrix generated using (1). The parameters of the ungapped Markov model are calculated using a so-called Gibbs sampler, written by Nielsen et al. [24 (link)].
The result of applying (1) is a prediction score for every residue of the query sequence. To reduce fluctuations, a smoothing window is applied to every position. It is made asymmetric in the N- and C- termini in order to conserve prediction examples.

5 different surface measures, calculated from the protein structure, were trained and tested for their ability to predict B-cell epitopes. These were variation of residue contact counts: Full sphere neighbor count (FS) [12] (link), Upper half-sphere neighbor count (UHS) [22] (link) and Half-sphere exposure as described in [22] (link) (HSE) and previously used for B-cell prediction in [18] (link). A residue were classified as neighbor to the query residue if the C_α - C_α distance were below k_sur. We furthermore tested the widely used relative surface accessibility (RSA) [21] (link) and a hybrid between neighbor count and RSA (Ta) by defining neighbor residues as residues holding any atom within T distance of any atom in the query residue. Scoring functions and parameters are listed in Table S3. Neighbor count in upper and lower half-spheres were calculated using the structural bio-python module developed by T. Hamelryck [32] (link), and surface accessibility calculated by DSSP using the standard 4 Å probe. The RSA were then obtained by dividing the surface accessibility with the maximum surface accessibility, calculated from the peptide GGXGG, where X is the amino acid in question. The optimal sphere radius k_sur, for the FS, UHS, RSA, and HSE, and the distance threshold T for Ta were estimated by a grid search using the grids; k_sur = {4,6…28 Å} and T = {4,6…28 Å}.

The filtered 4 proteins from the previous step were uploaded to SignalP- 5.0 server (Almagro Armenteros et al., 2019 (link)) to predict the location of signal peptides. Following that, the mature polypeptides were analyzed for their T and B cell epitopes through the Immune Epitope Database (IEDB) (Dhanda et al., 2019 (link)). Firstly, we mapped CTLs for the protein candidates using the HLA allele reference set, which provided more than 97% in terms of population coverage (Weiskopf et al., 2013 (link)), and the NetMHCpan EL 4.1 prediction tool (that was recommended by the IEDB database). Secondly, we mapped for HTLs against the HLA reference set to cover more than 99% in terms of population coverage (Greenbaum et al., 2011 (link)) and used IEDB recommended 2.22 as a prediction method. Furthermore, HTL peptides were assessed for their ability to induce several cytokines such as IFN-gamma (Dhanda et al., 2013b (link)), IL-4 (Dhanda et al., 2013a (link)), IL-10 (Nagpal et al., 2017 (link)), IL-6, and IL-13 (Jain et al., 2022 (link)). The last analysis for HTLs and CTLs was the conservancy prediction where multiple sequence alignment against the corresponding proteins in other reference sequences was employed to validate the conservancy of the selected epitopes. The last set of epitopes; namely BCLs were finally estimated through IEBD using the BepiPred-2.0 prediction method (Jespersen et al., 2017 (link)). Following prediction, the estimated epitopes were filtered based on the consideration of several characteristics such as the number of reacting alleles (to achieve a high population coverage percentage), conservancy percentage, and antigenicity score.

While continuous B cell epitopes, which were predicted in the above sections, are estimated through the primary amino acid sequence of a protein, another type of epitopes, which is conformational (or discontinuous) B cell epitopes, are predicted based on the 3D structure of the antigenic protein. For this purpose, the ElliPro Server (http://tools.iedb.org/ellipro) was utilized.

Sequence analysis included 80 nucleotide sequences of PfGARP from Thai isolates, one clinical isolate from Guinea (isolate MDCU32) and 18 publicly available complete gene sequences whose isolate names, country of origins and their GenBank accession numbers are as follows: 3D7 (Netherlands from West Africa, AL844501), CD01 (Congo, LR129686), Dd2 (Indochina, LR131290), FC27 (Papua New Guinea, J03998), FCC1/HN (Hainan in China, AF251290), GA01 (Gambia, LR131386), GB4 (Ghana, LR131402), KH1 (Cambodia, LR131418), KH2 (Cambodia, LR131306), HB3 (Honduras, LR131338), IGH-CR14 (India, GG6656811), IT (Brazil, LR131322), KE01 (Kenya, LR131354), ML01 (Mali, LR131481), SD01 (Sudan, LR131466), SN01 (Senegal, LR131434), TG01 (Togo, LR131450), and UGT5.1 (Vietnam, KE124372). Of these, the 3D7, FC27and FCC1/HN sequences were determined by Sanger dideoxy-chain termination method whereas the remaining isolates were assembled sequences from next-generation sequencing platforms (Supplemental Table S1). Sequence alignment was performed by using the CLUSTAL_X program, taken into account appropriate codon match in the coding region by manual adjustment to maintain the reading frame. The sequence from the FC27 strain was used as a reference^{6 (link)}. Searching for nucleotide repeats was performed by using the Tandem Repeats Finder version 4.0 program with the default option. Nucleotide diversity (π), the rate of synonymous substitutions per synonymous site (d_S) and the rate of nonsynonymous substitutions per nonsynonymous site (d_N) were determined from the average values of sequence differences in all pairwise comparison of each taxon and the standard error was computed from 1000 bootstrap pseudoreplicates implemented in the MEGA 6.0 program^{41 (link)}. Haplotype diversity and its sampling variance were computed by taking into account the presence of gaps in the aligned sequences using the DnaSP version 5.10 program^{42 (link)}. Natural selection on codon substitution was determined by using fast unconstrained Bayesian approximation (FUBAR) method in the Datamonkey Web-Server^{43 (link),44 (link)}. Neighbor-joining phylogenetic tree based on nucleotide sequences was constructed by using maximum composite likelihood parameter whereas maximum likelihood tree was built using Tamura-Nei model with the rate variation model allowed for some sites to be evolutionarily invariable. The Arlequin 3.5.2.2 software was deployed to determine genetic differentiation between populations, the fixation index (F_ST), using analysis of molecular variance approach (AMOVA) akin to the Weir and Cockerham’s method but taken into account the number of mutations between haplotypes^{45 (link)}. One hundred permutations were deployed to determine the significance levels of the fixation indices. Prediction of linear B cell epitopes in PfGARP was performed by using a sequence similarity to known experimentally verified epitopes from the Immune Epitope DataBase (IEDB) implemented in the BepiBlast Web Server^{11 (link)}. Furthermore, linear B cell epitopes were also predicted based on protein language models implemented in BepiPred-3.0^{12 (link)}. Potential HLA-class II-binding peptides were analyzed by using the IEDB recommended 2.22 algorithm with a default 12–18 amino acid residues option. The predicted HLA-class II-binding peptides were predicted based on the percentile rank < 10 and the IC₅₀ threshold for HLA binding affinity ≤ 1000 nM^{14 (link)}. The analysis mainly concerned the common HLA class II haplotypes among Thai populations with allele frequency > 0.1^{13 (link)}.