3.1. Structural Clustering
A total of 40, 22, and 38 crystal structures for BTK, JAK, and PARP were extracted. The crystal structures were clustered into four categories (Fig. 2). 6O8I (1.42 Å), 4RX5 (1.36 Å), 4RFZ (1.17 Å) and 5P9K (1.28 Å) were identified as the representative crystal structures of BTK with the highest resolution among each category (Additional file 2: Table S1), 6BBU (2.08 Å), 6N7A (1.33 Å), 6AAH (1.83 Å) and 4EHZ (2.17 Å) for JAK (Additional file 2: Table S2) and 4ZZZ (1.9 Å), 5WS1 (1.9 Å), 6I8T (2.1 Å) and 6NRH (1.5 Å) for PARP (Additional file 2: Table S3). The scatterplots of RMSD values for crystal structures versus residues were also calculated (Additional file 2: Figure S1).
3.2. Performance in Docking Calculations
We have previously studied that the compatible scoring functions by considering both reproducibility of ligand conformation and binding affinity prediction[34, 35]. The scoring function CHEMPLP was demonstrated to be appropriate in both docking power and scoring power and thus applied to scoring function (Additional file 2: Table S4).
The co-crystal ligand and prepared validation dataset was used to assess the docking power and discrimination power which presented by RMSD and P-value. The number of active compounds and decoys was listed in Table 1. Docking calculations were conducted by screening the validation dataset to discriminate the inhibitors from non-inhibitors. U-test (P-value) was used to evaluate the significance of the difference between the distributions of docking scores of the inhibitors and non-inhibitors.
The RMSD value below 2 Å is considered as a successful prediction of binding pose. Though, as shown in Table 2, in PARP the structure 6I8T is the unique that had succeed in pose prediction while four representative crystal structures of BTK and JAK have successfully reproduced the ligand conformation by re-docking the ligand as illustrated in Table 2 and Fig. 3. For BTK, each representative crystal structure also has The P-values suggested the significant difference of docking scores between inhibitors and decoys which means favorable discrimination power (Table 2, Fig. 4). The best P-values of three targets were 8.586E-74, 2.438E-67 and 6.260E-57 for BTK, JAK, and PARP, respectively. Apparently, the crystal structures of BTK performed best in distinguish the inhibitors and non-inhibitors. The docking performance of different target were different, even the different structures for the same target exhibited remarkable difference. On the basis of P-values, the representative complex 5P9K, 6N7A and 6I8T was chosen as the initial structure of BTK, JAK and PARP for MD simulations, respectively.
Table 1
Targets, actives and decoys used in this study
Target
|
BTK
|
JAK
|
PARP
|
actives
|
125
|
152
|
134
|
decoys
|
2500
|
3040
|
2680
|
EFmax
|
20
|
20
|
20
|
Table 2
The GOLD docking performance of representative crystal structures against BTK, JAK and PARP, respectively
Target
|
PDB ID
|
RMSD(Å)
|
P-Value
|
BTK
|
6O8I
|
1.6808
|
3.557E-71
|
4RX5
|
0.8877
|
9.094E-70
|
4RFZ
|
0.8824
|
2.954E-70
|
5P9K
|
1.6393
|
8.586E-74
|
JAK
|
6BBU
|
0.7786
|
1.846E-52
|
6N7A
|
0.3008
|
2.438E-67
|
6AAH
|
0.5392
|
1.093E-56
|
4EHZ
|
0.5708
|
3.234E-62
|
PARP
|
4ZZZ
|
2.0444
|
1.085E-32
|
5WS1
|
2.3612
|
1.079E-32
|
6I8T
|
0.9136
|
6.260E-57
|
6NRH
|
2.3873
|
2.046E-48
|
P-value is used to quantify the significance level of the two-tailed asymptotic significance of the Man-Whitney U test. The PDB ID of the bold font is selected as the initial structures of molecule dynamics simulations.
3.3. Molecular Dynamics Simulations
The crystal structures 5P9K, 6N7A and 6I8T were chosen as the initial structure to construct the MD simulations. For each complex, 40 ns MD simulations were enough for trajectory to achieve equilibrium. The structural deviation between the structures in the dynamic simulation process and the initial structure were calculated based on the backbone. As shown in Fig. 5, the RMSD of backbone for BTK, JAK and PARP is under 0.3 nm which means the system has equilibrated in the molecule dynamic simulation process. To tolerate structural differences in the stable phase of each protein, 20–40 ns trajectory was collected. The RMSD matrix between each two conformations was constructed, based on which can perform the subsequent cluster analysis. On this basis, the 20–40 ns trajectory was clustered, and the conformation with middle RMSD obtained in the cluster can be used as a representative structure (Fig. 6). For BTK, the structures 5P9K_3768, 5P9K_2526, 5P9K_2992 and 5P9K_2629 were recovered as representative MD snapshots. While the structures 6N7A_3398, 6N7A_2014, 6N7A_3726 and 6N7A_2010 as well as 6I8T_3041, 6I8T_3906, 6I8T_2072 and 6I8T_3274 were recovered as representative MD simulations for JAK and PARP respectively.
3.4. The Correlation Analysis of Docking Performance Based on Crystal Structures and Structures from MD Simulations
Docking performance was evaluated by docking scores which were obtained from molecular docking simulation. The docking scores of the known inhibitors and non-inhibitors of each target based on each representative crystal structure and structure from MD simulations were used to execute the correlation analysis including Pearson and Spearman ranking correlation analysis. When absolute value of the correlation coefficient (r or ρ) is greater than or equal to 0.8, it is considered that the two variables are highly correlated. While when less than 0.3 indicates weak correlation, basically irrelevant. The redundant structures whose correlation coefficient is high (༞0.9), which means the compared two structures possess highly similar docking results, were discarded. Then the Pearson correlation coefficients (r) and the Spearman ranking correlation coefficients (ρ) of the docking scores based on each two compared protein structures were calculated (Fig. 7). In detail, the results were concluded in Additional file 2: Table S5 and Additional file 2: Table S6.
3.5. Performance of Bayesian Models Based on Each Single Structure and Ensemble of Protein Structures
To constitute ensembles with different structures in size and element for exploring how different number and member of ensemble influence the virtual screening performance, the structures 6O8I and 5P9K_2629 in BTK, 6N7A_2681 in JAK as well as 6I8T_3906 in PARP are discarded based on the correlation analysis and 5-fold cross validation result for Bayesian Model. We have constructed the ensembles by following the principle: validation result is preferred, and the smaller-sized ensemble is the basis of the bigger. It needs to be noted that in consideration of the differentiation between crystal structures and MD simulations in validation result, the ensembles in which all members are crystal structures or MD simulations are validated. Based on these, nine panels are identified: Crystal Structure, MD Simulation, Two-size Ensemble, Three-size Ensemble, Four-size Ensemble, Five-size Ensemble, Six-size Ensemble, Seven-size Ensemble and Eight-size Ensemble. For all targets, the members of all ensembles are listed in Additional file 2: Table S7. Bayesian Models are created based on the docking scores of each protein structure and ensemble, and validated by leave-one-out cross-validation. The detailed 5-Fold Cross-Validation Results are summarized in Additional file 2: Table S8 and an external test set was used to test the performance of the model (Additional file 2: Table S9). The Bayesian score was extracted from the validation result and summarized in Table 3. For each single receptor, the virtual screening performance of Bayesian Model was quantified by the area under curve (AUC) of its receiver operating characteristic (ROC) plot. Furthermore, Four-size ensemble: all crystal structures combined and all MD simulations combined; the ensemble with maximal model score whose ROC curves are plotted as well.
It is presented in Table 3, with Bayesian score as the evaluation criterion that the protein structures of BTK performed well in both training set and external test with all ROC scores above 0.9. In turn, the protein structures of JAK and PARP exhibited worse performance with some ROC scores below 0.9 and even there are ROC scores below 0.8. In JAK the crystal structures outperformed the MD simulations while the reverse performance observed in BTK and PARP. The same trend was discovered in all crystal structures combined (Ensemble 7) and all MD simulations combined (Ensemble 8) ensembles in JAK and PARP. Comparing the single rigid structures with ensembles, universally, the ensembles regardless of different size give better Bayesian score which means that the Bayesian Models based on ensembles give more satisfactory predictions in identifying inhibitors and non-inhibitors. The best performed Bayesian Models in all targets are Ensemble 10 (Five-size Ensemble), Ensemble 10 (Five-size Ensemble) and Ensemble 9 (Four-size Ensemble) for BTK, JAK and PARP, respectively. Further, it is obviously that the Models’ performance does not always increase with the number of member. The best performing model usually happens to be a Four-size or Five-size Ensemble.
The ROC curves were used to illustrate the virtual screening performance of Bayesian Models (Fig. 8), the corresponding AUC were displayed in Table 4. It is consistent with the Bayesian score; the MD simulations outperformed the crystal structures in BTK and PARP. However, the crystal structures combined ensemble (Ensemble 7) was presenting more outstanding performance than the MD simulations combined ensemble (Ensemble 8). Lower correlation of docking scores between crystal structures compared with MD simulations may account for this. When different structures used for the docking calculations, the top ranking compounds are quite different. This elucidates the importance of choosing the appropriate protein structure in docking calculations, and then in ensemble constructing.
Table 3
Bayesian Score of the Bayesian Models Based on the Docking Scores of Each Single Representative Complex and Ensemble for BTK, JAK and PARP
Panel
|
BTK
|
JAK
|
PARP
|
Ensemble
|
ROC Score
|
Ensemble
|
ROC Score
|
Ensemble
|
ROC Score
|
Training Set
|
Test Set
|
Training Set
|
Test Set
|
Training Set
|
Test Set
|
Crystal
Structure
|
4RFZ
|
0.957
|
0.95
|
6BBU
|
0.862
|
0.856
|
4ZZZ
|
0.779
|
0.791
|
6O8I
|
0.96
|
0.948
|
6N7A
|
0.906
|
0.883
|
5WS1
|
0.775
|
0.828
|
4RX5
|
0.958
|
0.942
|
6AAH
|
0.87
|
0.854
|
6I8T
|
0.878
|
0.908
|
5P9K
|
0.961
|
0.97
|
4EHZ
|
0.89
|
0.87
|
6NRH
|
0.843
|
0.865
|
MD Simulation
|
5P9K_3768
|
0.951
|
0.948
|
6N7A_3398
|
0.828
|
0.783
|
6I8T_3041
|
0.879
|
0.906
|
5P9K_2526
|
0.967
|
0.952
|
6N7A_2014
|
0.864
|
0.796
|
6I8T_3906
|
0.853
|
0.92
|
5P9K_2992
|
0.968
|
0.983
|
6N7A_2681
|
0.796
|
0.768
|
6I8T_2072
|
0.841
|
0.87
|
5P9K_2629
|
0.963
|
0.941
|
6N7A_2010
|
0.855
|
0.834
|
6I8T_3274
|
0.867
|
0.934
|
Two-size Ensemble
|
Ensemble 1
|
0.975
|
0.973
|
Ensemble 1
|
0.913
|
0.894
|
Ensemble 1
|
0.905
|
0.925
|
Ensemble 2
|
0.974
|
0.986
|
Ensemble 2
|
0.884
|
0.836
|
Ensemble 2
|
0.892
|
0.934
|
Ensemble 3
|
0.978
|
0.992
|
Ensemble 3
|
0.912
|
0.887
|
Ensemble 3
|
0.912
|
0.938
|
Three-size Ensemble
|
Ensemble 4
|
0.979
|
0.982
|
Ensemble 4
|
0.916
|
0.896
|
Ensemble 4
|
0.897
|
0.923
|
Ensemble 5
|
0.976
|
0.985
|
Ensemble 5
|
0.879
|
0.835
|
Ensemble 5
|
0.895
|
0.938
|
Ensemble 6
|
0.981
|
0.991
|
Ensemble 6
|
0.917
|
0.885
|
Ensemble 6
|
0.918
|
0.939
|
Four-size Ensemble
|
Ensemble 7
|
0.981
|
0.98
|
Ensemble 7
|
0.926
|
0.911
|
Ensemble 7
|
0.893
|
0.92
|
Ensemble 8
|
0.977
|
0.982
|
Ensemble 8
|
0.874
|
0.829
|
Ensemble 8
|
0.895
|
0.933
|
Ensemble 9
|
0.981
|
0.981
|
Ensemble 9
|
0.919
|
0.899
|
Ensemble 9
|
0.919
|
0.948
|
Five-size Ensemble
|
Ensemble 10
|
0.983
|
0.988
|
Ensemble 10
|
0.926
|
0.912
|
Ensemble 10
|
0.916
|
0.949
|
Ensemble 11
|
0.982
|
0.989
|
Ensemble 11
|
0.924
|
0.9
|
Ensemble 11
|
0.914
|
0.942
|
Six-size Ensemble
|
Ensemble 12
|
0.983
|
0.987
|
Ensemble 12
|
0.924
|
0.901
|
Ensemble 12
|
0.914
|
0.946
|
Seven-size Ensemble
|
Ensemble 13
|
0.983
|
0.986
|
Ensemble 13
|
0.918
|
0.894
|
Ensemble 13
|
0.91
|
0.942
|
Eight-size Ensemble
|
Ensemble 14
|
0.982
|
0.985
|
Ensemble 14
|
0.914
|
0.888
|
Ensemble 14
|
0.906
|
0.939
|
Table 4
Virtual Screening Performance of Bayesian Models for BTK, JAK and PARP
Validation Result
|
BTK
|
JAK
|
PARP
|
Ensemble
|
AUC
|
Quality
|
Ensemble
|
AUC
|
Quality
|
Ensemble
|
AUC
|
Quality
|
4RFZ
|
0.916
|
Excellent
|
6BBU
|
0.827
|
Good
|
4ZZZ
|
0.751
|
Fair
|
6O8I
|
0.916
|
Excellent
|
6N7A
|
0.844
|
Good
|
5WS1
|
0.759
|
Fair
|
4RX5
|
0.91
|
Excellent
|
6AAH
|
0.836
|
Good
|
6I8T
|
0.859
|
Good
|
5P9K
|
0.93
|
Excellent
|
4EHZ
|
0.865
|
Good
|
6NRH
|
0.826
|
Good
|
5P9K_3768
|
0.923
|
Excellent
|
6N7A_3398
|
0.794
|
Fair
|
6I8T_3041
|
0.86
|
Good
|
5P9K_2526
|
0.924
|
Excellent
|
6N7A_2014
|
0.805
|
Good
|
6I8T_3906
|
0.83
|
Good
|
5P9K_2992
|
0.931
|
Excellent
|
6N7A_2681
|
0.749
|
Fair
|
6I8T_2072
|
0.816
|
Good
|
5P9K_2629
|
0.92
|
Excellent
|
6N7A_2010
|
0.817
|
Good
|
6I8T_3274
|
0.854
|
Good
|
Ensemble 7
|
0.974
|
Excellent
|
Ensemble 7
|
0.919
|
Excellent
|
Ensemble 7
|
0.905
|
Excellent
|
Ensemble 8
|
0.972
|
Excellent
|
Ensemble 8
|
0.858
|
Good
|
Ensemble 8
|
0.901
|
Excellent
|
Ensemble 10
|
0.979
|
Excellent
|
Ensemble 10
|
0.919
|
Excellent
|
Ensemble 9
|
0.922
|
Excellent
|