Due to the scarcity of suitable antibody-antigen pairs, computational docking was used to generate some of the data in the training and testing data set. The ClusPro (Boston University) [6-9] and Rosetta (RosettaCommons) [10-12] web servers were used to create a data set of paired antibody-antigen complexes for machine learning. Both ClusPro and Rosetta were used for protein-protein molecular docking. Rosetta uses the SnugDock (RosettaCommons) algorithm [10]. The Swiss-PdbViewer (Swiss Institute of Bioinformatics) [13] was used to examine the resulting protein complex structures.

A total of 50 antibody-antigen complexes were selected randomly from the Protein Data Bank (PDB) [14]. The antibody-antigen complexes were separated by using a Perl script to produce PDB-formatted files as well as sequences for antibodies and antigens. CDRs were located by using the Rosetta antibody modeling web server. Antigens were docked with a range of antibodies by using ClusPro (used only to determine orientation), followed by Rosetta’s antibody docking program, SnugDock. In order to keep computation times manageable, not all antibodies were docked. Instead, 10 to 14 antibodies were randomly selected to be docked with each antigen in order to find the best orientation. The resulting complexes were submitted to the Rosetta SnugDock web server in order to calculate the best interface score. This produced structures for between 10 and 14 complexes per antigen, which, when added together with the original antibody-antigen complex, totaled 11 to 15 complexes per antigen. Altogether, 50 antigens were docked with 600 antibodies. An example of a resulting complex is shown in Figure 1.

The Rosetta interface scores were used as estimates of binding affinity in order to identify cognate antibody-antigen pairs to be used as input for machine learning. Complexes with interface scores of higher than −8.0 were classified outright as complexes with poor binding, and those with interface scores of lower than −9.0 were classified outright as complexes with good binding. For complexes with scores that ranged between −8.0 and −9.0, the docking clusters and positions were examined visually by using SwissDock (Swiss Institute of Bioinformatics). If the top 10 models had their antibodies and antigens in similar relative positions and the structures showed sensible interaction patterns, the pairs were classified as having a good binding affinity.

Rosetta interface scores have been used previously as classifiers to determine binding affinity based on docking results (eg, in an antibody-antigen cross reactivity study [15]).

Additional data were extracted from the Coronavirus Antibody Database (CoV-AbDab) [16]—a database of antibodies against coronaviruses, including SARS-CoV-2, SARS-CoV-1, and MERS-CoV (Middle East respiratory syndrome–related coronavirus). Data (2674 rows) were extracted from the CoV-AbDab on February 14, 2021. After filtering out incomplete data, 2031 rows remained, with each row corresponding to an antibody. The information extracted comprised the antibody names, their binding antigens, and their heavy and light variable region sequences, including the locations of the third CDRs (CDR3s). Each of the variable region sequences were searched against the international ImMunoGeneTics information system database [17] in order to identify the locations of the first CDRs (CDR1s) and second CDRs (CDR2s) from the heavy and light chains. Since a row may contain information about an antibody's interactions with multiple antigens, the data were further split into multiple rows, with each row containing information about the interaction between 1 antibody and 1 antigen.

Additional features were calculated for the sequences, as follows. The isoelectric point for each CDR was calculated by using the Bachem peptide calculator analysis tool (Bachem Holding AG) [18]. The average hydrophilicity of each CDR was also calculated by using the Bachem peptide calculator.

B cell epitopes were predicted by using the IEDB (Immune Epitope Database) antibody epitope prediction analysis tool [19].

The resulting data set can be downloaded from GitHub [20] and is structured with the following column headings: H chain CDR1 sequence, H chain CDR2 sequence, H chain CDR3 sequence, L chain CDR1 sequence, L chain CDR2 sequence, L chain CDR3 sequence, Hydrophilicity of L CDR1, pI of L CDR1, Hydrophilicity of L CDR2, pI of L CDR2, Hydrophilicity of L CDR3, pI of L CDR3, Hydrophilicity of H CDR1, pI of H CDR1, Hydrophilicity of H CDR2, pI of H CDR2, Hydrophilicity of H CDR3, pI of H CDR3, Antigen Epitope, Rosetta Docking score, Antigen, and Docking result.

A weighted K-nearest neighbor (K-NN) classification algorithm [21] for predicting antibody-antigen binding affinity was implemented in Python. The program can be downloaded from GitHub [20].

For each antigen, the 11 to 15 antibodies that were docked were labeled as “good affinity” or “low affinity,” on the basis of the docking results. Machine learning was then performed, using the sequences of both antigens and antibodies.

Neighbors were determined by using the string distances between the CDR1, CDR2, and CDR3 amino acid sequences of different antibodies. Weights were calculated from distances, so that nearer neighbors were considered to have more weight, as detailed below.

For every antigen, the class (good affinity or low affinity) was learned by using the K-NN method, using a training subset (N − 1) of the labeled antigen-antibody sequence pairs and using the CDR string distances as features. The model performance was then evaluated on the remaining antigen-antibody sequence pair that was not used for training (leave-one-out cross-validation).

In order to ensure that the K-NN pairs only included pairs with the same antigen, a fixed penalty of 1000 was added to the distances between antibody-antigen pairs involving different antigens.

The similarity between antibodies was measured via a comparison of their CDRs. Each antibody has a heavy chain and a light chain, and each chain contains 3 CDRs. The distance between 2 antibodies was calculated as the Euclidean distance between their CDR distance vectors, as shown in the following equation (equation 1):

where (q_i – p_i) represents the string distance between the CDR_i of antibody q and the CDR_i of antibody p.

The Python code is given in Multimedia Appendix 1.

Two different CDR distance calculation methods were tested and compared; one was based on sequence identity, and the other used the BLOSUM62 matrix, as detailed below.

For the identity-based distance measure, pairs of equivalent CDRs were compared with each other based on their Levenshtein string distances [22], as shown in the following equation (equation 2):

Cost=0 for a_i=b_i, Cost=1 for a_i≠b_i

The Levenshtein distance only accounts for amino acid identity when it is used for comparing sequences. A more biologically significant distance measure needs to take into account the different properties of amino acids, which means that some amino acid substitutions are more likely to be accepted in an interaction than others. The BLOSUM62 substitution matrix [23] was used as a proxy for amino acid similarity in the Levenshtein distance calculation. Although the BLOSUM matrices were designed to reflect evolutionary conservation, they can provide an estimate of similarity in interaction potential [24].

The Levenshtein distance was calculated as per equation 2, using the following cost function:

For a_i=b_i, Cost=0

where S_ij, S_ii, and S_jj are obtained from the BLOSUM62 matrix.

The following columns from the data set were used to train the model for leave-one-out cross-validation: H chain CDR1 sequence, H chain CDR2 sequence, H chain CDR3 sequence, L chain CDR1 sequence, L chain CDR2 sequence, L chain CDR3 sequence, Antigen, and Docking result. The trained model was then evaluated on its ability to predict the docking results from the other columns.

A random forest machine learning algorithm incorporating the previous K-NN results was also used for predicting antibody-antigen binding classification. The isoelectric point and net charge at neutral pH (7.0) for each CDR were used as additional features, in addition to the BLOSUM62-derived CDR distances, for training the random forest. Binding was predicted by combining the votes from each of the features, and each individual feature contributed 1 vote, according to the nearest neighbor predictions based on each feature.

The following columns from the data set were used for training the random forest: String distance (calculate by KNN method), Hydrophilicity of L CDR1, pI of L CDR1, Hydrophilicity of L CDR2, pI of L CDR2, Hydrophilicity of L CDR3, pI of L CDR3, Hydrophilicity of H CDR1, pI of H CDR1, Hydrophilicity of H CDR2, pI of H CDR2, Hydrophilicity of H CDR3, pI of H CDR3, Antigen, and Docking result. The trained model was then evaluated on its ability to predict the docking results from the other columns.

Each feature was considered as an individual decision tree and contributed 1 vote. For example, the isoelectric point of the CDR1 of an antibody’s heavy chain was considered as 1 feature, and the K-NN method was used, as previously described, to find the results of this decision tree. Altogether, there were 13 decision trees, and each tree used the K-NN method to determine its vote, for a total of 13 votes. The final decision was determined based on a simple majority vote. The best results were obtained when the whole forest (all 13 decision trees) took part in the vote.

The performance of the K-NN and random forest learners was evaluated by using leave-one-out cross-validation on an antigen basis. For each of the 57 antigens, a training data set was constructed by removing 1 row, that is, 1 antibody-antigen pair, from the data set. After training with the remaining antibodies that bound to this antigen, model performance was evaluated based on the removed antibody. The process was repeated until all 5041 antibody-antigen pairs were tested. Model accuracy was calculated as the ratio of the number of correctly predicted antibody-antigen pairs over the total number of pairs in the data set.

A total of 600 antibody-antigen complexes were generated via the computational docking of 50 antibody structures with 50 antigen structures. In addition, a total of 4441 antibody-antigen pairs were extracted from the Cov-AbDab. The composition of this section of the data set is shown in Table 1.

In total, the data set contained 5041 antibody-antigen pairs comprising 1157 antibodies and 57 antigens.

The antigen-antibody binding classification methods were evaluated by using leave-one-out cross-validation. For a K value of 2 nearest neighbors, the K-NN method, when the Levenshtein distance was calculated based on sequence identity, achieved an accuracy of 81%. A slight improvement (accuracy of 82%) was observed when using the BLOSUM62 matrix to calculate the Levenshtein string distance.

Different K values were also evaluated when the Levenshtein distance was calculated based on the BLOSUM62 matrix. A K value of 2 provided the best accuracy. For a K value of 1 nearest neighbor, the accuracy was 80%. For a K value of 3, classification accuracy dropped to 79%.

For the random forest predictions, votes were used as the classification prediction results. The accuracy was highest when the whole forest was considered, in which case each feature contributed to the classification results. The performance of the random forest method was best (accuracy of 80%) when all 13 features—the Levenshtein string distance and the isoelectric point and net charge at neutral pH (7.0) for each CDR—took part in the final votes.