Deep Learning of High-Order Interactions for Protein Interface Prediction

A protein complex is composed of two proteins, known as the ligand protein and the receptor protein. Suppose the ligand protein has $n_{l}$ amino acids and the receptor protein has $n_{r}$ amino acids, there are $n_{l}\times n_{r}$ possible amino acids pairs. Protein interface prediction problem aims at predicting if there exists an interaction within each amino acid pair. Following the existing work (Fout et al., 2017), we represent proteins as graphs, where nodes represent amino acids and edges indicate affinities between amino acids. Formally, we define the feature matrix of the ligand protein as $\hat{H}^{l}=[\hat{\mathbf{h}}^{l}_{1},\hat{\mathbf{h}}^{l}_{2},\cdots,\hat{\mathbf{h}}^{l}_{n_{l}}]\in\mathbb{R}^{d\times n_{l}}$ and the feature matrix of the receptor protein as $\hat{H}^{r}=[\hat{\mathbf{h}}^{r}_{1},\hat{\mathbf{h}}^{r}_{2},\cdots,\hat{\mathbf{h}}^{r}_{n_{r}}]\in\mathbb{R}^{d\times n_{r}}$ where $d$ denotes that each node has a $d$-dimensional feature vector.

The existing work (Fout et al., 2017) formulates it as a binary classification problem that predicts the interaction between each node pair separately. Specifically, for the the $i$-th node in $\hat{H}^{l}$ and the $j$-th node in $\hat{H}^{r}$, it concatenates the corresponding feature vectors $\hat{\mathbf{h}}^{l}_{i}$ and $\hat{\mathbf{h}}^{r}_{j}$ , then uses dense layers as a binary classifier to determine whether the two nodes interact with each other. However, such a formulation ignores high-order pairwise interactions that the interaction prediction of an amino acid pair may also depend on the information of other amino acid pairs. In addition, converting amino acids to graphs loses the sequential information of the original amino acid chains.

To address these issues, we incorporate the sequential information and formulate the protein interface prediction as a 2D dense prediction problem. First, given the node feature matrices $\hat{H}^{l}$ and $\hat{H}^{r}$, we propose the sequential modeling (SM) to restore the sequential information, which results in the order-preserved node feature matrix $H^{l}=[\mathbf{h}^{l}_{1},\mathbf{h}^{l}_{2},\cdots,\mathbf{h}^{l}_{n_{l}}]\in\mathbb{R}^{d\times n_{l}}$ for the ligand protein and $H^{r}=[\mathbf{h}^{r}_{1},\mathbf{h}^{r}_{2},\cdots,\mathbf{h}^{r}_{n_{r}}]\in\mathbb{R}^{d\times n_{r}}$ for the receptor protein. Next, we propose the high-order pairwise interaction (HOPI) to generate a 3D tensor $\mathcal{Q}\in\mathbb{R}^{n_{l}\times n_{r}\times c}$ where each $\mathcal{Q}_{ij}\in\mathbb{R}^{c}$ denotes the feature combination of the $i$-th node in $H^{l}$ and the $j$-th node in $H^{r}$. Finally, based on the tensor $Q$, the protein interface prediction problem predicts a 2D matrix $O\in\{0,1\}^{n_{l}\times n_{r}}$. Each element of $O$ can be either 0 or 1. For location $O_{i,j}$, 1 means there exists an interaction between $i$-th amino acid of the ligand protein and $j$-th amino acid of the receptor protein, while 0 indicates there is no interaction. Since the predictions $O$ is generated based on the whole tensor $Q$, both sequential information and high-order pairwise interactions are incorporated.

Both structural and sequential information are important for studying properties of protein complexes. As shown in Figure 1, representing proteins as graphs can well-convey the structural information. It is a popular way since we can employ graph neural networks to pass, transform and aggregate structural information across graphs (Ying et al., 2018). However, the sequential information from the original amino acid chains is lost when converting proteins to graphs because of order-invariant property of graphs. Such a sequential structure is the primary structure of a protein, which is unique to other proteins and defines important functions of the protein.

To overcome this limitation, we propose the sequential modeling to preserve the sequential information of the original input amino acid chains for a given protein. Formally, given an input protein with $m$ amino acids, we first record the original sequential order set $I=(1,2,\cdots,m)$. Then we formulate it as graphs and map each node in the graph with the index in $I$. Next, we employ graph neural networks to learn node features, denoted as $\hat{H}=[\hat{\mathbf{h}}_{1},\hat{\mathbf{h}}_{2},\cdots,\hat{\mathbf{h}}_{m}]\in\mathbb{R}^{d\times m}$, where $d$ is the dimension of a node feature vector. Then we reorder the feature matrix $\hat{H}$ based on the order set $I$ as

Then the node feature vectors in the new feature matrix $H$ has a consistent order with the original amino acid sequential order. In this way, feature matrix $H$ successfully captures both structural information from the protein graph and the sequential information from the amino acid chains. We believe such a reordering operation helps capture complex inherent relationships among amino acids, thereby resulting in more accurate interface prediction.

Protein interface prediction aims at determining interactions between two amino acids from different proteins. The protein structure is the 3D arrangements in amino acid chains and always folds into specific spatial conformations to enable biological functions (Pauling et al., 1951). It is possible that any two amino acids from different proteins can interact with each other. Existing methods (Fout et al., 2017; Sanchez-Garcia et al., 2019; Townshend et al., 2019) predict the interaction for each amino acid pair separately. One amino acid is picked from the ligand protein graph and the other is from the receptor protein graph. The features of the two amino acids are concatenated and passed to dense layers for binary classification. However, high-order context for amino acid pairs are ignored, which can help extract the important high-level interaction patterns. Hence, we propose the high-order pairwise interactions to learn complex interaction patterns for interface prediction.

Suppose we have the sequence-preserved node feature matrix $H^{l}=[\mathbf{h}^{l}_{1},\mathbf{h}^{l}_{2},\cdots,\mathbf{h}^{l}_{n_{l}}]\in\mathbb{R}^{d\times n_{l}}$ for the ligand protein and $H^{r}=[\mathbf{h}^{r}_{1},\mathbf{h}^{r}_{2},\cdots,\mathbf{h}^{r}_{n_{r}}]\in\mathbb{R}^{d\times n_{r}}$ for the receptor protein. We compute a third-order tensor $\mathcal{Q}\in\mathbb{R}^{n_{l}\times n_{r}\times c}$. Each $\mathcal{Q}_{ij}\in\mathbb{R}^{c}$ in $\mathcal{Q}$ is the transformation of $\mathbf{h}^{l}_{i}$ and $\mathbf{h}^{r}_{j}$. It can be computed by either summation of $\mathbf{h}^{l}_{i}$ and $\mathbf{h}^{r}_{j}$ or the concatenation of $\mathbf{h}^{l}_{i}$ and $\mathbf{h}^{r}_{j}$. The proposed HOPI allows the tensor $\mathcal{Q}$ to store structural information, sequential information, and inherent high-order pairwise interactions information. Then we employ convolutional neural networks (CNNs) to perform 2D dense predictions based on the tensor $\mathcal{Q}$. Stacking several CNN layers extracts high-level pairwise interaction patterns from a region containing a subsequence from the ligand protein, a subsequence from the receptor protein, and inherent high-order pairwise interactions from the two subsequences. Finally, the output $O\in\{0,1\}^{n_{l}\times n_{r}}$ indicates the interactions between any possible amino acid pairs. Note that the prediction of $O_{i,j}$ depends not only on $\mathcal{Q}_{ij}$ but also on all feature interactions within its receptive field.

The overall architecture is illustrated in Figure 2. Given a ligand protein and a receptor protein, GNNs are used to pass, transform and aggregate the structural information in protein graphs. All GNN layers are shared by the two proteins. Then the proposed sequential modeling performs reordering to preserves the sequential information of the original amino acid chains for both proteins. Next, the proposed high-order pairwise interaction method produce a 3D tensor containing feature interactions for all amino acid pairs. Finally, the tensor is passed to 2D CNNs for 2D dense prediction. The output map contains interaction predictions for all amino acid pairs. Note that any modern CNN architecture such as ResNet (He et al., 2016), UNet (Ronneberger et al., 2015) or DeepLab (Chen et al., 2017) can be flexibly integrated into our framework to perform dense prediction and the whole system can be trained end-to-end.

We employ graph neural networks to aggregate structural information. Suppose a node $\mathbf{h}_{i}$ in the protein graph has $n$ nodes in its neighborhood. The neighboring node feature matrix is $H_{i}\in\mathbb{R}^{d^{N}\times n}$, and the neighboring edge feature matrix is $E_{i}\in\mathbb{R}^{d^{E}\times n}$, where $d^{N}$ is the dimension of node feature vectors and $d^{E}$ is the dimension of edge feature vectors. We first aggregate both node and edge features from neighborhood as

where $W^{N}\in\mathbb{R}^{d^{N}\times d^{N}}$, $W^{E}\in\mathbb{R}^{d^{N}\times d^{E}}$, and $\mbox{tanh}(\cdot)$ is an element-wise operation. $W^{N}$ and $W^{E}$ are used to perform linear transformation on the neighboring node features and edge features, respectively. The node feature matrix $H_{i}$ and edge feature matrix $E_{i}$ can be treated as a set of node vectors $[\mathbf{h}_{i1},\cdots,\mathbf{h}_{ij},\cdots,\mathbf{h}_{in}]$ and a set of edge vectors $[\mathbf{e}_{i1},\cdots,\mathbf{e}_{ij},\cdots,\mathbf{e}_{in}]$, respectively. Note that node vectors in $H_{i}$ and edge vectors in $E_{i}$ are in a consistent order. An edge $\mathbf{e}_{ij}$ links the center node $\mathbf{h}_{i}$ to the corresponding node $\mathbf{h}_{ij}$. The achieved matrix $Z_{i}$ aggregates information from both nodes and edges within the neighborhood. Each vector $\mathbf{z}_{ij}\in\mathbb{R}^{d^{N}}$ contains information from the node $\mathbf{h}_{ij}$ and the edge $\mathbf{e}_{ij}$.

We use two methods to transform neighborhood information to the center node. The first one is to simply perform average on all neighboring vectors $\mathbf{z}_{ij},j=1,\cdots,n$, namely neighborhood average (NeiA). Another approach is neighborhood weighted average (NeiWA), which essentially assigns relatively larger weights to these important vectors while smaller weights to the ones that are less important, then performs the weighted average on the neighboring nodes and edges. We introduce the two approaches below.

Protein interface prediction is to determine whether there are interactions between amino acids from two different proteins. Essentially, it investigates cross-protein pairwise interactions. The interactions can be partly determined by features and inherent properties of amino acids in both proteins. We name an interacted amino acid pair as a positive sample and a non-interacted pair as a negative sample. Generally, the number of positive samples is significant less than that of the negative samples in a protein complex. Hence, it causes the data unbalance problem. To address this issue, we propose to incorporate in-protein interaction information to increase the number of positive examples, and hence improve the predictions of cross-protein interactions.

Specifically, we propose to use HOPI to capture both in-protein and cross-protein pairwise interactions. Given the sequence-reserved node feature matrix $H^{l}\in\mathbb{R}^{d\times n_{l}}$ for the ligand protein and $H^{r}\in\mathbb{R}^{d\times n_{r}}$ for the receptor protein, the achieved tensor $\mathcal{Q}$ is expanded to the size of $(n_{l}+n_{r})\times(n_{l}+n_{r})\times c$. An illustration of the tensor $\mathcal{Q}$ is provided in Figure 3. Either of the two regions $n_{l}\times n_{l}$ and $n_{r}\times n_{r}$ contains in-protein structural and sequential information, and in-protein pairwise interactions. The two regions $n_{l}\times n_{r}$ are same as those introduced in Section 3.3, which contain cross-protein interactions. 2D CNNs are performed such that in-protein structural and sequential information, in-protein pairwise interactions, and cross-protein pairwise interactions are all captured for interface prediction.

We define a training sample as a pair of amino acids and is labeled by their interaction. All samples in a protein complex can be treated as a sub-epoch. During one iteration, the network is trained on a part of samples in a protein complex. In this way, two subgraphs on the ligand and the receptor graphs are updated when aggregating neighboring structural information using GNNs. A small patch is updated on the tensor $\mathcal{Q}$ when performing dense prediction using CNNs. A sub-epoch is finished when all training samples in one graph complex are used. And an epoch is finished when all training samples in all graph complexes in the dataset are passed to the network.

When considering the cross-protein pairwise interactions only, a common situation is that the number of positive samples is much less than that of the negative samples in a protein complex. Down-sampling on the negative samples is usually required to reach a reasonable predefined positive-negative (PN) ratio, which allows the training of the network. However, all negative samples are kept in the prediction phase. One way for data augmentation is to incorporate in-protein pairwise interactions. By doing this, some in-protein positive samples could be added in the training process. This is expected to compensate the small PN ratio in the cross-protein pairwise interactions.

We use three datasets to evaluate our proposed methods. They all come from the popular Docking Benchmarks, which include several protein-protein docking and binding affinity benchmarks. The first one is generated from Docking Benchmarks version 5 (DB5) (Vreven et al., 2015; Fout et al., 2017). It is the largest and the most recent dataset which contains complicated protein structures for protein interface prediction. The dataset was originally split into training, validation and test sets (Fout et al., 2017; Townshend et al., 2019; Sanchez-Garcia et al., 2019). The training set contains 140 complexes, the validation set contains 35, and the test set contains 55 complexes. Each complex is a pair of ligand and receptor proteins. A data sample is a pair of amino acids from different proteins and their interaction. The total number of positive samples in the dataset is $20857$, which is much less than the number of negative samples. The PN ratio is around 1:1000. The Docking Benchmark version 4 (DB4) (Hwang et al., 2010) and Docking Benchmark version 3 (DB3) (Hwang et al., 2008) are two subsets of the DB5. The DB4 contains 175 complexes and 16004 positive samples, and the DB3 contains 127 complexes and 12335 positive samples in total.

All dataset are mined from the Protein Data Back (Berman et al., 2000). Each protein has original sequential information from amino acid chains. The numbers of amino acid in protein sequences ranges from tens to thousands. The protein structures are obtained from X-ray crystallography or biological mutagenesis experiments (Afsar Minhas et al., 2014). In protein graphs, nodes are amino acid and edges are affinities between nodes. both nodes and edges contain features from protein structures and sequences.

We use the same node features as in (Fout et al., 2017; Afsar Minhas et al., 2014). The node features are computed based on different properties of amino acid. The residue depth is defined as the minimal distance for an amino acid to the protein’s surface. It’s normalized in the range from 0 to 1 and has been demonstrated to carry valuable information for amino acid interactions. The amino acid composition defines the count of a specific amino acid in the direction and opposite direction of the side chain for the amino acid of interest. The threshold along two directions is the minimal atomic distance of 8A. The amino acid composition varies dramatically among amino acids, which is vital to determine the properties of an amino acid. The protrusion index for an amino acid is a collection of statistics of protrusion values for all atoms along its side chain. These features deliver important inherent structural information and properties for the amino acid of interest. They are combined and concatenated together in a consistent order, which results in the total number of node features to be 76. We use the same edge features as in the work (Fout et al., 2017). Each edge feature vector contains 2 features. One is the normalized distance between two amino acids and the other is the angle between two normal vectors for the two amino acid planes.

Recently a larger dataset Database of Interacting Protein Structures (DIPS) is created by (Townshend et al., 2019). An amino acid is represented as 4D grid-like data, which contains spatial information at the atom level and types of all atoms in the amino acid. However, the structural information is not considered in the dataset, thus the proteins can not be represented as graphs.

The baseline methods could be grouped into three categories, these are, the state-of-the-art conventional machine learning method BIPSPI (Sanchez-Garcia et al., 2019), the CNN-based method SASNet(Townshend et al., 2019) and the GNN-based methods DCNN (Atwood and Towsley, 2016), NGF (Duvenaud et al., 2015), DTNN (Schütt et al., 2017) and NEA (Fout et al., 2017). Particularly, the GNN-based baselines use different graph neural architectures for node feature learning, but use the same dense layers as binary classifiers to predict the interaction for each pair of amino acid separately.
BIPSPI is the abbreviation for xgBoost Interface Prediction of Specific-Partner Interactions. The method combines both structure and sequence features and uses Extreme Gradient Boosting (Chen and Guestrin, 2016) with a novel scoring function for protein interface prediction.
SASNet is the Siamese Atomic Surfacelet Network, which uses only spatial coordinates and types of all atoms in amino acids and voxelizes all amino acids into a 4D-grid manner. The first three dimensions deliver the spatial information of an amino acid and the last dimension is the one-hot representation of types for all atoms in the amino acid. The paired two amino acid representations are then passed to 3D CNN with weights sharing, followed by concatenation operation and dense layers for binary classification to decide whether the two amino acids interact with each other.
DCNNs is diffusion-convolutional neural networks for graph-structured data applying diffusion-convolution operators $k$ times ($k$-hops) for node feature learning. A diffusion-convolution operator scans a diffusion process for each node. For a node of interest, $k$ diffusion-convolution operators gather information from all nodes that each of those can connect to the node of interest through $k$ steps. Then several dense layers are used as a binary classifier to predict the interactions of two nodes.
NGF is the commonly used graph convolutional networks, which first aggregates node information in neighborhood by multiplying the adjacency matrix to the node feature matrix, and then performs linear transformations on node features, followed by a nonlinear function for node feature learning.
DTNN is deep tensor neural networks, which aggregates both node and edge information in neighborhood. Linear transformations are applied to node features and edge features separately. After that, element-wise multiplication on the corresponding node features and edge features is performed to achieve the final feature vector for the node of interest. Intuitively, edges serve as gates to help control information from the corresponding nodes.
NEA is Node and Edge Average, the state-of-the-art GNN-based method on the used DB5 dataset. Similar as DTNN, it performs aggregation and linear transformation on both nodes and edge features in neighborhood. Then node and edge features are averaged and summed together, followed by a residual connection and nonlinear activation to generate node features.

We use the same data splitting as in all the baseline methods (Fout et al., 2017; Townshend et al., 2019; Sanchez-Garcia et al., 2019) for the DB5 datasets. For the DB4 and DB3 datasets, we first randomly split each dataset with a ratio of 6:2:2 for the training, validation and test samples, then fix the splitting in all experiments. We first only consider the cross-protein pairwise interactions for training. Similar to the work (Fout et al., 2017), we keep all the positive examples and perform down-sampling on the negative samples, resulting in the PN ratio of 1:10 during the training phase. We maintain the original PN ratio in the validation and test phases.

For three GNN-based baselines, different numbers of GNN layers are designed and explored in (Fout et al., 2017). We also conduct experiments using the same numbers of GNN layers for fair comparisons. Our proposed GNN layer has two variants, neighborhood average and neighborhood weighted average. We conduct experiments on both for clear comparisons. For high-order pairwise interaction, we perform concatenation on the feature vectors of two nodes from different proteins. Several residual blocks are employed for dense prediction. A residual block contains two 2D convolutional layers, the first of which is followed by ReLU as the activation function. The number of intermediate channels of the first convolutional layer is set as a hyperparameter. The identity map of the input is summed to the output of the second convolutional layer, followed by ReLU to generate the final output.

We use the grid search to tune hyperparameters. The search space for all hyperparameters is provided in Table 1. Adam Optimizer (Kinga and Adam, 2015) is employed for training and ReLU is used as the activation function. Each experimental setting is conducted to run 10 times with different random seeds. All hyperparameters are tuned based on the validation set. Optimal hyperparameters are tuned on one run and used across all the 10 runs.

As positive examples and negative examples are not balanced, we use Receiver operating characteristic (ROC) curve for evaluation. Specifically, we calculate Area Under the ROC Curve (AUC) based on the ROC curve for each complex. Then median AUC (MedAUC) for all the complexes in the test sets is used to evaluate the performance of different models. The used MedAUC can grantee very large or very small proteins will not have dramatic effects on the performance on the whole dataset.