We calculate the sequence similarity between each pair of input peptides in three ways. The input for each method is the list of peptide sequences, and the output is a 2D matrix containing the pairwise similarity scores. The first sequence similarity score is calculated using a BLOSUM matrix, where each pair of amino acids is assigned an integer value based on the relative frequencies of amino acids (18). The BLOSUM62 matrix values are calculated based on amino acid sequence alignments with less than 62% identity. In our method, the similarity between two peptides is defined as the sum of the BLOSUM62 values for the amino acid pair at each position of the peptides, which is a common method of calculating sequence similarity. Secondly, we calculated the pairwise similarity between peptides using the similarity matrix calculated by HLAthena (19), which calculates the entropy at each peptide position based on the entire dataset and uses the entropy to weight the importance of each peptide position when calculating similarity based on the PMBEC similarity matrix. Lastly, we calculated the Hamming distance between two peptides, defined as the number of amino acid positions that differ between the two peptides (i.e., AAAA and AAAB have a Hamming distance of 1 because only position 4 differs). The combination of these three similarity metrics was empirically observed to give the best results.

The pHLA pairwise similarity is calculated starting from a dataset of pHLA structures. The output is another 2D matrix of pairwise similarity scores. Depending on the source of the structures, their reference frame may be different, so we first align the structures using the align function of PyMOL (20). Then, we extract the solvent-accessible surface value from each structure using the program MSMS with a density of 3.0 and a probe size of 1.5Å (21). MSMS computes the surfaces as a triangular mesh, which we then downsample to a resolution of 1.0Åusing pymesh (22). We then annotate the vertices of this mesh with biochemical features, specifically the electrostatic potential, hydrophobicity, and hydrogen bond potential. The electrostatic potential of the surface is calculated using APBS (23). The hydrophobicity of each amino acid in the pHLA is assigned based on the Kyte-Doolittle scale (24), and assigned to each surface point based on the closest amino acid. Finally, the hydrogen bond potential at each point is calculated based on the free hydrogens of the closest amino acid residues. We calculate the hydrogen bond potential using the data preparation method of MaSIF (25), based on an orientation-dependent hydrogen bonding potential (26). In brief, the hydrogen bond potential at a vertex is calculated based on the vertex’s distance and angle from potential hydrogen donors (polar hydrogens) and potential acceptors (nitrogen or oxygen). The potential ranges between -1 (hydrogen bond acceptor) and +1 (hydrogen bond donor).

To account for the T-cell only interacting with a specific part of the pHLA complex, we define the TCR-interacting region as a round patch centered on the peptide bound to the HLA. The center of the patch is calculated by finding the closest surface point to the peptide’s center of mass. A circular patch is extracted by selecting the vertices in the neighborhood of the center point up to 16 edges away from the center, as defined by the triangular mesh. The surface patch mesh is converted to a point cloud where each point is a vertex from the mesh, and each point is annotated with the biochemical features.

To calculate the similarity between pHLAs, we first perform a pairwise alignment of the point clouds using the Iterative Closest Point (ICP) algorithm (27). This alignment uses only geometric information and does not take into account the annotated biochemical features. The ICP is an iterative procedure that aligns a source point cloud S to a target cloud T in 3D space by iterating over three main steps. The first is to create a corresponding point set C by matching points in the source point cloud to points on the target point cloud within a distance of ϵ = 2 Å. By taking points within distance ϵ of each other, we account for the fact that there may be an unequal number of points in point clouds S and T . The p s i l o n ball approach is a standard procedure in ICP (27). The second step is to calculate the rotation and translation that will best minimize the distance D between each corresponding point pair (i.e., to find the best transformation to align each source point to its corresponding target point). For a pair of point clouds S and T and the set of corresponding points C , the distance D is defined as

Where S i is the vector of x , y , and z coordinates of point i of point cloud S . The third step of ICP is to transform the source points using the rotation and translation found in the previous step. These three steps are repeated, recreating the corresponding point set for each iteration until convergence: when the change in D is less than 1.0 e − 06 , or until 30 iterations have been performed. ICP results in a final distance D and the indices of the corresponding points for each point cloud.

We define a new score function to calculate the similarity between aligned point clouds using the geometric coordinates and the biochemical features. We expand the ICP distance D so that each vertex in the point clouds has six dimensions, to include not only the geometric coordinates x , y , and z , but also the biochemical features electrostatic potential, hydrogen bond potential, and hydrophobicity. We also account for the size of the corresponding point set, as a low number of corresponding points indicates that the source and target point clouds are not well aligned. Our distance score D 2 is defined as

This score is calculated between each pair of point clouds, resulting in a 2D matrix where is element s , t is the score between point clouds s and t .

Each similarity calculation method described above (i.e., BLOSUM62, HLAthena, Hamming distance, structural and biochemical) produces a 2D matrix, where element i , j is the similarity between peptides/pHLAs i and j . Each matrix is normalized by subtracting the mean of the matrix and dividing it by the standard deviation. Then, all matrices are summed element-wise, producing a 2D matrix. This matrix provides the similarity score between each pair of pHLAs.

The pairwise similarity scores produced by our method can be used to cluster the peptides. To validate the method, we performed agglomerative clustering using scikit-learn (28). We performed clustering with ward linkage or average linkage. We used agglomerative clustering to create 2 clusters (with parameters n_clusters=2 and distance_threshold=None) to represent the two possible clusters “cross-reactive peptides” and “non-cross-reactive peptides”. We also ran agglomerative clustering to create any number of clusters (with n_clusters=None and distance_threshold=0) and used the result to build a dendrogram to visualize the distance between the different peptides.

We also validated the method using K-nearest-neighbors (KNN) clustering (29). KNN is supervised clustering, meaning the true label of each peptide is known, except the peptide we are attempting to label based on its nearest neighbors. We used k values between 1 and 8 and used KNN to label each peptide based on all the other peptides in the dataset.

We tested our method on five datasets, as explained below. The full list of peptides in each dataset is provided in the Supplementary Data.

The PepSim scoring method is available through a web server interface (see Figure 1 ). Users may upload a dataset in PDB format, including the peptide that they wish to target. The target peptide will be used as the reference for the final ranking of peptides based on similarity. After submission and execution, users can visualize the peptides in a dendrogram based on agglomerative clustering of the peptides based on the similarity scores. The peptides are also visualized in a 2D scatter plot created by using non-metric multidimensional scaling (NMDS) (36). The users also receive a ranked list of the peptides based on their similarity to the given target. If they want to perform offline analysis, including clustering, users can download all the results, which include the 2D array of pairwise similarity scores.

The PepSim web server is implemented using Docker (37), with the backend implemented in Django (38). The submitted jobs are managed by a distributed task queue by Celery (39). The webserver is currently hosted on a virtual machine in the Owl Research Infrastructure Open Nebula (ORION) VM Pool on Rice University Campus.

To test the accuracy of our similarity score, we performed the same pipeline on five datasets. On dataset 1, containing 61 peptides that are recognized by the same T-cell and 119 negative decoys, we performed agglomerative clustering with the similarity scores to create two clusters. The resulting clustering had a sensitivity of 98.36% and a specificity of 96.64%. A visualization of the clustering results, as well as the true clustering, can be seen in Figures 2A, B . This clustering produced 1 false negative (ASDPMNHYY), and 4 peptides that have not been experimentally determined to produce a T-cell response (ELDPTNMTY, DSDPTGTAY, ELDPDNETY, ELDPNNAVY).

We also performed agglomerative clustering on dataset 3, containing 28 HCV peptides and 45 other viral peptides. One of the HCV peptides causes no T-cell response and the other 27 cause some level of T-cell response. As seen in Figures 3A, B , all 28 HCV peptides are clustered in the same cluster, and the decoys are in the other cluster. Given that all but one HCV peptide trigger a T-cell response, this clustering produces only one false positive (G3-18). This exemplifies that our similarity score can differentiate between peptides from different origins, as all the decoys are different viral peptides.

We performed a similar experiment on datasets 2, 4, and 5. In these cases, we did not specify the number of clusters and produced dendrograms for a full visual of the clusters. Dataset 2 contains the same 28 HCV-derived peptides as dataset 3. We see in Figure 4A that the final clustering separates the peptides with a high or intermediate T-cell response from most of the peptides with low or no response. The peptides with high or intermediate responses are split into different clusters, but within their separate clusters, they are clustered together. The dendrogram defines two clusters, but if we split the larger cluster based on the dendrogram branches within it, we completely separate the set of high-response peptides from the low and no response groups. The other set of high and intermediate response peptides are in the other cluster in the dendrogram, with one pHLA that produces a low T-cell response (G1-03).

Datasets 4 and 5 are smaller datasets with a low number of cross-reactive peptides. Dataset 4 shows a complete separation of cross-reactive and non-cross-reactive peptides, as seen in Figure 5 . The cross-reactive peptides are very similar in sequence. Dataset 5 contains two cross-reactive peptides: HEV-1527 and MYH9-478. As seen in Figure 6 , HEV-1527 and MYH9-478 are not separated into a different cluster from the negative decoys, but they are on the same branch of the dendrogram. If you select the target peptide HEV-1527, then the closest peptide is MYH9-478, and vice versa. Given that this method would start with a target peptide, then we can accurately identify the cross-reactive peptide from this dataset.

In 2015, Mendes used multivariate statistical methods to perform structure-based prediction of T-cell cross-reactivity among the 28 viral peptides in dataset 2 (8). Using their method, they clustered the peptides into two distinct clusters. The peptides that trigger a high T-cell response are placed in one cluster and the peptides with a low response are in the other cluster. One of the peptides with an intermediate response in the cluster with the high response peptides, and the other two are with the low response peptides. In contrast, with our method, the peptides with a high response are split into two different clusters, and the peptides with an intermediate response are clustered with the peptides with a high response.

Dataset 3, being an expansion of dataset 2, was previously studied in 2011 by Antunes et al. (7) Only 10 of the HCV peptides were included in the dataset (A0201_0031, A0201_0051, A0201_0052, A0201_0053, A0201_0054, A0201_0055, A0201_0056, A0201_0057, and A0201_0058), along with 45 other peptides. All 10 of the HCV peptides were clustered together, along with five other peptides not derived from HCV (A0201_0014, A0201_0083, A0201_0076, A0201_0095, and A0201_0073). In contrast, our method results in all the HCV peptides clustered together with only one false positive, A0201_0033, which is an HCV peptide that causes no T-cell response. We can also look at our results on dataset 2, where peptides A0201_0051-58 are clustered together with no other clusters, and A0201_0031 is in the other cluster.

To examine the effectiveness of the structural and biochemical analysis, we repeated the experiments with two variations of the similarity score. We defined one variation of the score to use only sequence analysis, and the other variation to use only the structural and biochemical analysis.

In dataset 1, when we use only sequence analysis, the final agglomerative clustering (using Ward linkage) has a sensitivity of 96.72% (2 false negatives) and specificity of 97.48% (3 false positives), as seen in Figures 2C, D . Compared to the similarity with the structural and biochemical analysis included, there is one more false negative and one less false positive when we remove the structural and biochemical analysis.

When we use the similarity score calculated from only the structural and biochemical analysis, the sensitivity increases to 100%, but the specificity decreases to 44%, as seen in Figures 2E, F .

We also used K-nearest-neighbors to determine the effectiveness of the similarity score. As seen in Table 1 , we achieved the highest sensitivity and specificity when we use all the components of the similarity score (sensitivity of 100% and specificity of 98.3%). We also achieve high accuracy with the variations of the similarity score, including using only the structure and biochemical analysis

In dataset 3, removing the structural and biochemical analysis from the similarity score does not change the final clustering when using Ward or Average linkage, as seen in Figures 3C, D . Similarly, using only the structural and biochemical analysis results in one additional false positive when using Ward linkage (A0201_0014). Figures 3E, F shows the 2D embedding of the peptides based on the structural and biochemical similarity score and compared to Figures 3A, B the clusters are not as separated.

We also used K-nearest-neighbors to determine the effectiveness of the similarity score. As seen in Table 2 , we achieved the highest sensitivity and specificity when we use all the components of the similarity score (sensitivity of 100% and specificity of 95.6%). We also achieve high accuracy with the variations of the similarity score, including using only the structure and biochemical analysis.

justification=centering

Dataset 2 provides more interesting results. When the structural and biochemical analysis is removed from the similarity score, the dendrogram clustering shows that the high and intermediate responders are still separated into two different clusters, as seen in Figure 4B . One of the clusters has three false positives compared to the one false positive in Figure 4A .

When using only structural and biochemical analysis, there is one cluster of all peptides with a high or intermediate T-cell response, and the other peptides with a high response are separated from each other ( Figure 4C ). Therefore, the clusters are less accurate when only using the structural and biochemical analysis when compared to the complete score. However, the peptides derived from different viral genotypes are clustered together.

We also used K-nearest-neighbors to determine the effectiveness of the similarity score. As seen in Table 3 , we achieved the highest sensitivity and specificity when we use all the components of the similarity score (sensitivity of 100% and specificity of 92.9%). We also achieve high accuracy with the variations of the similarity score, including using only the structure and biochemical analysis.

Datasets 4 and 5 show little change when removing the structural and biochemical analysis and when removing the sequence analysis. In dataset 4, the four cross-reactive peptides and four non-cross-reactive peptides are separated into different clusters with 100% accuracy. In dataset 5, again HEV-1527 and MYH9-478 are not separated into a different cluster from the negative decoys, but they are on the same branch of the dendrogram.