PASSer: prediction of allosteric sites server

The data used in the current work was collected from the Allosteric Database (ASD) [26]. There are a total of 1946 entries information of allosteric sites with different proteins and modulators. To ensure data quality, 90 proteins were selected using previous rules [11]: protein structures with either resolution below 3Å or missing residues in the allosteric sites were removed; redundant proteins in the rest of the data that have more than 30% sequence identity were filtered out. The names and IDs of the 90 proteins extracted from the PDB [27, 28] are listed in table S1 (is available online at stacks.iop.org/MLST/2/035015/mmedia).

FPocket algorithm [29] is used to detect pockets from the surface of the selected proteins. A pocket is labeled as either 1 (positive) if it contains at least one residue identified as binding to allosteric modulators or 0 (negative) if it does not contain such residues. Therefore, a single protein structure may have more than one positive label. A total of 2246 pockets were detected with 119 pockets being labeled as allosteric sites. There are 19 features calculated by the FPocket as shown in table S2.

Pearson correlation coefficient (PCC) measures the linear correlation of two variables. Given a pair of variables X and Y as $\{(x_1, y_1), \ldots, (x_n, y_n)\}$, PCC (r_X,Y ) is calculated as

$$\begin{equation}r_{X, Y} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i=1}^n (y_i - \bar{y})^2}}\end{equation} \tag{ 1 }$$

where n is the sample size and $\bar{x},\kern1.5pt \bar{y}$ are sample means. PCC has a value between −1 and +1. Positive value represents a positive correlation, and negative value represents a negative correlation. The absolute value indicates the degree of correlation. The larger the absolute value, the stronger the correlation.

XGBoost is an ensemble learning method that combines several decision trees in sequence.

Let $D = \{(x_i, y_i) (|D| = n, x_i \in R^m, y_i \in R^n) \}$ represents a dataset with m features and n labels. The jth decision tree in XGBoost predicts a sample $(x_i, y_i)$ by

$$\begin{equation} g_j(x_i) = w_q(x_i) \end{equation} \tag{ 2 }$$

where w_q is the leaf weights of this decision tree. The final prediction of XGBoost is given by the summation of predictions from each decision tree:

$$\begin{equation} \hat{y_i} = \sum_{j = 1}^M g_j(x_i) \end{equation} \tag{ 3 }$$

where M is the total number of decision trees. To overcome overfitting introduced by decision trees, the objective function in XGBoost is composed of a loss function l and a regularization term Ω:

$$\begin{equation} \text{obj} (\theta) = \sum_{i = 1}^N l(y_i, \hat{y_i}) + \sum_{j = 1}^M \Omega(f_i) \end{equation} \tag{ 4 }$$

where $\Omega(f) = \gamma T + \frac{\lambda}{2}\sum_{l = 1}^T \omega_l^2$ with T represents the number of leaves and γ, λ are regularization parameters.

During training, XGBoost iteratively adds new decision trees. The prediction of the tth iteration is expressed as

$$\begin{equation} \hat{y_i}^{(t)} = \hat{y_i}^{(t-1)} + g_t(x_i). \end{equation} \tag{ 5 }$$

Correspondingly, the objective function of the tth iteration is

$$\begin{equation} \text{obj}^{(t)} = \sum_{i = 1}^N l(y_i, \hat{y_i}^{(t-1)} + g_t(x_i)) + \Omega(f_i). \end{equation} \tag{ 6 }$$

XGBoost introduces both first derivative and second derivative of the loss function. By applying Taylor expansion on the objective function at second order, the objective function of the tth iteration can be expressed as

$$\begin{align} \text{obj}^{(t)} &\simeq \sum_{i = 1}^N [l(y_i, \hat{y_i}^{(t-1)}) + \partial_{\hat{y}^{(t-1)}} l(y_i, \hat{y}^{(t-1)}) f_t(x_i) \nonumber\\ & \quad+ \frac{1}{2} \partial_{\hat{y}^{(t-1)}}^2 l(y_i, \hat{y}^{(t-1)}) f_t^2(x_i)] + \Omega(f_i). \end{align} \tag{ 7 }$$

XGBoost can predict the labels of sample data with corresponding probabilities. For one pocket, XGBoost outputs the probability of this pocket being an allosteric pocket. This pocket is labeled as positive (allosteric) if the predicted probability is over 50% or negative otherwise.

There are only 5.3% of positive labels in this binary classification job, which means that the input data is highly imbalanced. To focus more on the limited positive labels, XGBoost uses 'scale_pos_weight', a parameter for controlling the balance of positive and negative weights. A typical value of this weight equals to the number of negative samples versus the number of positive samples.

In the current study, the maximum tree depth for base learners and the weights for positive labels were fine-tuned while keeping other parameters as default values. The XGBoost algorithm is implemented using Scikit-learn package [30] version 0.23.2.

The GCNN in this work follows this formula [22]:

$$\begin{equation} H^{(l + 1)} = \text{ReLU}(\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}H^{(l)}W^{(l)}) \end{equation} \tag{ 8 }$$

where H^(l) and H^(l + 1) represents the lth and l + 1th layer, respectively. $H^{(l)} \in \mathbb{R}^{N \times D}$, where N is the number of nodes and D is the number of features. Rectified linear unit (ReLU(x) = max(0, x)) is used as the activation function. W^(l) denotes the weight matrix in the lth layer. D and A represent degree matrix and adjacent matrix, respectively, with $\tilde{D}_{ii} = \sum_j \tilde{A}_{ij}$. Renormalization (indicated by the ∼ symbol) is applied for the undirected graph G where each node is added with a self-connection. Therefore, $\tilde{A} = A + I_N$ where I_N is the identity matrix.

A graph readout is calculated through the average of node features for each graph:

$$\begin{equation} h_g = \frac{1}{|\mathcal{V}|} \sum_{v \in \mathcal{V}} h_v \end{equation} \tag{ 9 }$$

where h_g is the readout result of graph g and h_v is the node feature in node v. $\mathcal{V}$ represents the nodes in graph g.

An example of 1-layer GCNN model is shown in figure 1. A graph is first fed into a convolution layer. The in-degree and out-degree of a node refer to the number of edges coming into and going out from that node, respectively. In-degree of each node was calculated as the node feature. Graph feature is calculated as the average of node features in the readout layer with ReLU activation function. The output was further fed into a linear classification layer g, which predicts the probability of being an allosteric pocket. Previous research [31] has shown the limitations of 1-layer GCNN. In the current study, atomic graphs of each protein pocket are constructed and fed into a 2-layer GCNN model. This model consists of two graph convolution layers, each with 256 dimensions of a hidden node feature vector, followed by a readout layer and a linear classification layer. The node degree is used as the initial node feature. The in-degree is the same as the out-degree in an undirected atomic graph. Graph representation is calculated as the average of node representations. The linear classification layer outputs the probabilities of pockets being allosteric sites.

Zoom In Zoom Out Reset image size

To overcome the potential limitation in training GCNN with an imbalanced dataset, the ratio between negative labels and positive labels was fine-tuned. Specifically, in each protein, allosteric pockets (positive samples) were fully used, while non-allosteric pockets (negative samples) were partially used. Each negative sample was randomly selected, and the total number of non-allosteric pockets equals the number of allosteric pockets times the ratio value.

In constructing atomic graphs, the threshold of bond distance was also fine-tuned. Each atom was considered as a node in the atomic graph and the pairwise distances between atoms were calculated. If the distance is below a specified threshold, an edge is constructed connecting the two related nodes. Thus, the distance threshold controls the degree of local connectivity.

The GCNN model is implemented using deep graph library package [32] version 0.4.3.

For binary classification, the results can be classified as table 1. Various indicators were used to quantify model performance: precision (TP / (TP + FP)) measures how well the model can predict real positive labels; accuracy ((TP + TN) / (TP + FP + FN + TN)) measures the overall classification accuracy; recall (or sensitivity, TP / (TP + FN)) and specificity (TN / (TN + FP)) together measure the ability to classify TP and TN; F1 score (2 * precision * recall / (precision + recall)) is the weighted average of precision and recall. The higher the values of these indicators, the better the model's performance.

A receiver operating characteristic (ROC) curve was applied as another indicator to test model performance in binary classification. ROC curve is plotted as TP rate against FP rate with different threshold settings. The area under curve (AUC) is calculated for quantification. The upper limit for AUC value is 1.0. A dummy model should have an AUC value of 0.5.

The ranking information based on the predicted probabilities of protein pockets was also considered as an important indicator. Specifically, for each allosteric pocket, the ranking (predicted probability of being an allosteric site in descending order) was recorded and categorized as the first, second, third or other positions. The ranking result indicates how likely the allosteric pockets can appear in the top positions among all detected pockets in the same protein. A model with good performance should rank an allosteric pocket in the top positions.

The distribution of 19 features is shown in figure S1. While some exhibited long-tail distributions such as features 1 (Score) and 2 (Druggability score), data normalization is unnecessary since XGBoost does not require normal distribution. Instead, XGBoost, like other tree-based models, only focuses on the order, and whether features are normalized or not does not affect the prediction results. Violin plots are shown in figure S2 to better distinguish the feature distribution of allosteric sites and non-allosteric sites.

The correlation matrix between these features is shown in figure S3. Several features exhibited high correlations. For example, features 3 (Number of alpha spheres), 4 (Total SASA), 5 (Polar SASA), 6 (Apolar SASA), and 7 (Volume) are highly correlated with each other. Features 17 (Alpha sphere density) and 18 (Center of mass) are also strongly correlated with these five features.

For each feature value in a pocket, the inner ranking refers to the ranking position of this feature among values of other pockets in the same protein. Similar to the definition of inner ranking, the overall ranking refers to the ranking position of this feature value among values of other pockets in the overall dataset. Both sets of ranking features were normalized. The correlations between the original features and these two ranking feature sets were calculated and shown in figure S3. The high negative values indicate strongly negative correlations between the original features and the ranking features. While feature rankings were calculated and applied as additional features in a previous study [13], in the current dataset, the high correlation indicated that the ranking features provided little additional information and thus were discarded.

XGBoost model can overcome the limitation of data imbalance by controlling the weight difference between negative labels and positive labels. This parameter was fine-tuned along with the maximum depth of trees. The results are plotted in figure 2. Two sets of parameters reached high F1 scores, and both were selected in the final model. The final XGBoost model is composed of two models, each with one set of parameters. The results in any given pocket are the averaged results predicted by these two models.

Zoom In Zoom Out Reset image size

The results of the fine-tuned XGBoost model are listed in table 2. Compared with the reference results, XGBoost model exhibited higher accuracy, precision, specificity, and ROC AUC values with comparable results in recall and F1 scores. Therefore, XGBoost model performs well in allosteric site prediction.

Unlike XGBoost, GCNN models suffer from imbalanced dataset. To address this problem, the ratio between negative labels and positive labels was evaluated first. The results are plotted in figure 3(A). A ratio of 2 (number of negative labels: number of positive labels = 2:1) was selected. The distance threshold was further fine-tuned, and the results are plotted in figure 3(B). $10\ {\mathring{\mathrm A}}$ was selected as the distance cutoff when constructing atomic graphs.

Zoom In Zoom Out Reset image size

The results of the fine-tuned GCNN model with 10 independent runs are listed in table 2. Compared with XGBoost, GCNNs are less effective in classifying allosteric sites. However, it is expected that combining XGBoost and GCNN will result in better performance than either model.

The ensemble learning model is composed of both XGBoost model and GCNN model. For a given pocket, physical properties are calculated and fed into the XGBoost model; a representative atomic graph is fed into the GCNN model. The final result is calculated as the averaged probability of these two models. This final model contains both the physical properties and topological features of protein pockets. The combined results are listed in table 2. Compared with the XGBoost model, model ensembling leads to a 6.00% increase in recall, a 0.82% decrease in precision, and a 2.89% increase in F1 score. The AUC ROC value also had a 1.89% increase.

For each protein, the identified pockets are ranked based on the predicted probabilities. Overall, 60.7% of allosteric pockets are predicted as the first position, while 81.6% among the top 2 and 84.9% among the top 3. In other words, if a pocket is an allosteric pocket, there is a probability of 84.9% that it can be predicted in the top 3 among all detected pockets in the same protein. The prediction results, together with values reported in other studies, are listed in table 3. It should be noted that the types and amounts of proteins in the test set are different from each other.

To test this ensemble learning method, two proteins not in the dataset (table S1) were used. The predicted allosteric pockets of these two proteins are illustrated in figure 4. These two proteins represent two different types of allosteric proteins. The second PDZ domain (PDZ2) is a dynamics-driven protein in human PTP1E protein which undergoes allosteric process upon binding with peptides [33]. The light-oxygen-voltage domain of Phaeodactylum tricornutum Aureochrome 1a (AuLOV) is a conformational-driven allosteric protein [34]. AuLOV is a monomer in the dark state and undergoes dimerization upon blue light perturbation [35]. In both cases, our prediction model ranks the allosteric sites as the top 1 with probabilities of 45.14% and 89.46%, respectively. This indicates that this model is capable of predicting both dynamics-driven and conformational-driven allosteric proteins. Apparently, the probability for the dynamics-driven allosteric protein is much smaller than the one of the conformational-driven allosteric protein, which is not unexpected.

Zoom In Zoom Out Reset image size

A web server based on the allosteric prediction method developed in this study is implemented using a Python web framework, Django. JSmol [36] is a JavaScript implementation of the Jmol package and is embedded in the web page for protein and pocket visualization. Web pages are rendered using Bootstrap. This server is named as Protein Allosteric Sites Server (PASSer). A workflow of PASSer is outlined in figure 5.

Zoom In Zoom Out Reset image size

An example of input and output of PASSer is displayed in figure 6. Users can submit a PDB ID if available or upload a custom PDB file as shown in figure 6(A). By default, all chains in the protein are analyzed. Prediction results are displayed as two parts: top 3 pockets with the highest probability rendered with the protein structure (figure 6(C)) and their probabilities (figure 6(B)). For each pocket, the corresponding residues can be retrieved by clicking the 'Show Residues' texts. Protein structure is visualized using JSmol. Each pocket is either displayed upon clicking its 'Load pocket' icon or hidden by clicking its 'Hide pocket' or overall 'Reset' icons.

Zoom In Zoom Out Reset image size

A command line interface (CLI) is provided to facilitate potential developments. Similar with the web usage, this CLI can take either a PDB ID or a local PDB file for predictions.