Identifying DTIs can be regarded as a supervised prediction task to predict whether a pair of counterparts interact with each other or not in the drug-target networks. In this study, the benchmark dataset was taken from (He et al., 2010). There are mainly two reasons, 1) The information about the DTIs was collected from the DrugBanks, BRENDA, SuperTarget, and KEGG BRITE databases, which included four main drug target proteins of G Protein-coupled receptors (GPCR), enzymes (Ezy), ion channels (Chl), and nuclear receptors (NR). 2) In recent years, many researchers have been proposed to predict DTIs, which are based on this benchmark dataset, and hence will facilitate the comparison under the same condition. It can be summarized as follows: { S = S G P C R − D r u g + S C h l − D r u g + S E z y − D r u g + S N R − D r u g S G P C R − D r u g = S G P C R − D r u g + ( 630 ) + S G P C R − D r u g − ( 1240 ) S C h l − D r u g = S C h l − D r u g + ( 1372 ) + S C h l − D r u g − ( 2744 )               S E z y − D r u g = S E z y − D r u g + ( 2719 ) + S E z y − D r u g − ( 5438 )         S N R − D r u g = S N R − D r u g + ( 82 ) + S N R − D r u g − ( 164 )                                 (1)

There are 4,803 drug-target pairs in positive subsets, 2,719 for enzymes, 1,372 for ion channels, 630 for GPCRs, and 82 for nuclear receptors. Negative samples are randomly synthesized by separating each target and drug in S⁺, and none of them appear in the corresponding positive dataset. The proportion of positive samples and negative samples was set as 1:2. For comparison with previously published papers, both our positive and negative samples are consistent with He et al. (He et al., 2010)

Check390 is a dataset constructed by Hu et al. It contains 130 pairs of positive samples from the KEGG database, and 260 negative samples generated using the above method (Hu et al., 2016). Each pair in Check390 cannot be found in S .

In this article, we construct a novel model for DTIs based on large-scale pre-trained Bidirectional Encoder Representations from Transformers (BERT) and the fully connected neural network-based module called the BRL block. Figure 1 showers an overview of the DTIs model. The model has four modules: feature engineering, feature extraction, optimization, and decision-making. Firstly, in the feature engineering module, we use the auto-encoder ProtBert model, which is pre-trained on data from UniRef100 containing 216M protein sequences, to generate embedding vectors for protein sequences. As a result, the proteins can be represented via 1024-D vectors (dimensionality of the features extracted by the ProtBert model). Drug molecular fingerprints are represented by 128-D vectors through semi decomposition process discrete wavelet transform (DWT). Secondly, the 1152-D vectors (a concatenation of protein sequence feature and drug feature) are fed into the feature extraction model to generate interaction information through the first BRL block and CNN Afterwarderwards, in the decision-making module, the second BRL block is used to map interaction features into a unified vector space. The optimization module contains a contrastive loss and a cross-entropy loss. The contrastive loss is used to calculate the interaction information (generated by CNN block), which can reduce the distance between samples with the same label, and increase the distance between samples with different labels, while the cross-entropy loss is computed as the loss of second BRL block, bathes are used to adapt weights in the module during the learning process by minimizing the total loss. At the end of model, we can obtain the interaction score (generated by a softmax layer after second BRL block, and range from 0-1), the pair is interaction if the prediction score is   > 0.5 .

The CNN block includes a convolution layer, a rectified linear unit activation (ReLU), and a max-pooling layer. Instead of using multi-channels, we applied one channel only (Peng et al., 2018). In the convolution layer, apply a convolution kernel with a window size of h*k to extract the DTIs features, then use the rectified linear unit activation function and performed max-pooling to get the most useful interaction feature from the feature matrix subsequently. Through this block, an output of input xis formulated as:   v = m a x p o o l ( f ( w ⋅ x + b ) ) (12) where w ∈ R h k , which is applied to a window of h = 18 , k = 64 to produce a new feature; b ∈ R is a bias term and f is a non-linear function.

The BRL is built as a special block in the neural network, where data is normalized and then mapped into a specific vector space. This block consists of three layers: a batch-norm layer (BN), a leaky rectified linear activation layer (Leaky ReLU), and a linear layer (Pedregosa et al., 2011).

The input data x is first Batch-normalized, which serves to increase the learning rates further, remove the dropout layer, and apply other modifications afforded by the batch normalization (Ioffe and Szegedy, 2015); then input to the Leaky ReLU activation layer, and finally linearly mapped. BRL block can mathematically be represented as: X = L i n e a r   ( L e a k y R e L U ( B N ( x ) ) ) = W × ( L e a k y R e L U ( B N ( x ) ) ) + B (13) where x is the input data, the BN transform is applied independently to each dimension of x , W is the weight of the linear layer, and B is the bias of the linear layer. The first BRL block and CNN block are used for capturing both global and local information to represent the drug-protein pair; the second BRL block is used for predicting DTIs.

The BRL block was implemented with PyTorch (version 1.6.0), and a fully connected layer was used for the linear mapping. The parameters of the first BRL block were set as: the number of input neurons and the batch normalized dimensions dimension were both 1,152, and the number of output neurons was set to 128. The parameters of the second BRL block were set as 192 (128-D from the first BRL block and 64-D from the CNN block), and two respectively. A softmax layer is applied after the second BRL block, which is used to generate the prediction score. Other hyperparameters used default values in Pytorch. The source code for the related methods is available on a GitHub repository at: https://github.com/Jane4747/DTI-BERT.

In this frame, given two vectors v 1 and v 2 , input them into the same network in turn, the network will map the inputs to the new vector space where the similarity between two inputs can be evaluated by the distance measure function. Here, Euclidean distance was served as the distance measure, denoted as D ( v 1 ,   v 2 ) : D ( v 1 , v 2 ) = ‖ v 1 − v 22 ‖ 2 (14)

To make the samples of the same interaction label closer, and the distance between different labels as far as possible, the contrastive loss was applied as the loss function of the CNN network:   ℒ 1 ( v 1 , v 2 , Y ) =   1 2 ( 1 − Y ) D ( v 1 , v 2 ) 2 + 1 2 Y { m a x ( 0 , m − D ( v 1 , v 2 ) ) 2 } (15) where Y = 0 if sequences v 1 and v 2 , have the same label and Y = 1 if they are different, m > 0 is a margin. In other words, the margin defines a radius, and dissimilar pairs contribute to the loss function only if their distance is within the radius.

In this study, the second BRL block was used to convert the representation vector v to binary category outputs, the backpropagation algorithm was used to update network parameters, and the cross-entropy loss function was selected as the loss function of the second BRL block: ℒ 2 ( v , Y ) =   − Y   l o g ( D ( v ) ) − ( 1 − Y ) l o g ( 1 − D ( v ) ) (16)

Therefore, the loss function of the DTI-BERT model is: ℒ ( v 1 , v 2 , Y , Y 1 , Y 2 ) = ℒ 1 ( v 1 , v 2 , Y ) + ℒ 2 ( v 1 , Y 1 ) + ℒ 2 ( v 2 , Y 2 ) (17) where Y 1 and Y 2 are the labels of v 1 and v 2 .

We implemented our model using Python three and Pytorch (version 1.6.0). Optimizer, training epochs and batch size are set with “Adam”, 70 and 64, respectively. In our work, the optimizing function, “Adam”, use its default parameters value. All codes and trained models can be found via https://github.com/Jane4747/DTI-BERT.

The determination of a pair belongs to an interactive drug-target pair or non-interactive drug-target pair, is in the case of single-label classification. The metrics such as accuracy (ACC), sensitivity (Sn), Specificity (Sp), strength (str, the average of Sn and Sp) and Matthew’s correlation coefficient (MCC) are frequently used. The specific formulas are as follows: { A c c = T P + T N T P + T N + F P + F N S n =   T P T P + F N S p =   T P T P + F P S t r =   S p + S n 2 M C C =   T P × T N − F P × F N ( T P + F P ) ( T P + F N ) + ( T N + F P ) ( T N + F N ) (18) where TP represents the true positive, FN the false negative, TN the true negative, FP the false positive.

On the protein representation task, auto-encoder models (word2vec and BERT) with different model parameters scales were tested. For the drug representation task, a variety of algorithms in various fields, including natural language processing (word2vec), graph (node2vec and GCN), and signal processing (DWT) were tested.

We evaluated the BERT_Mean + DWT feature extraction method and compared it with several other classic protein and drug feature extraction methods, such as Pr ord2vec (a 64-D vector is obtained to represent the protein, it was extracted by an un-supervised word2vec model and implicated important biophysical and biochemical information (Yang et al., 2018; Zhang et al., 2020), BERT_First (the first row of P B E R T is obtained to represented protein, it is a 1024-D vector) (Nambiar et al., 2020), FP2_Word2vec (Jaeger et al., 2018), drug_Node2vec (Grover and Leskovec, 2016; Tetko et al., 2020), drug_Word2vec (Zhang et al., 2020; Zheng et al., 2021), drug_GCN (Chen et al., 2020). Figures 3–6 show the Matthews correlation coefficient (MCC) for the datasets S G P C R − D r u g , S C h l − D r u g , S E z y − D r u g , and S N R − D r u g obtained for each approach in CNN + BRL classifier via 10-fold cross validation.

It was found that BERT_Mean for the proteins and DWT for drugs can improve the performance of the classifier greatly in four datasets. The BERT_Mean + DWT increased capacity for identifying DTIs compared to the using BERT_First, PRO_Word2vec, drug_Node2vec, drug_Word2vec, and drug_GCN, and BERT_Mean can find the most compact and informative features subsets which are deeply hidden in protein sequences. It is showed that word2vec for protein sequences and GCN for drugs in DTIs tasks, could also obtain good prediction results on three datasets ( S G P C R − D r u g , S E z y − D r u g , and S N R − D r u g ), which inspires us that different protein representation methods need to consider different drug molecule representation methods, which need to be determined experimentally.

In order to test the performance of the BRL + CNN and compare it with the existing machine learning methods, we use the same benchmark dataset (listed in Eq. 1) and the same BERT_Mean + DWT feature as the input of the prediction model. The proposed BRL + CNN predictor and other commonly used classifiers provided by the Scikit-learn library, like Multi-Layer Perceptron (MLP) with two hidden layers (Pedregosa et al., 2011) and gradient boosting tree-based ensemble method called LightGBM (LGB) (Ke et al., 2017), were tested via 10-fold cross-validation, the results are listed in Table 1. It was found that the proposed BRL + CNN predictor in this article has better performance than other classifiers in all metrics.

To further demonstrate the power of the DTI-BERT predictor, we compared it with some existing methods. There are some new models for identifying DTIs trained with the datasets established by He et al. (He et al., 2010). For example, Hu et al. proposed a deep learning-based method to predict DTIs by using the information of drug structures and proteins sequences (Hu et al., 2019), this CnnDIT predictor has better prediction performance in predicting DTIs, and it has its own web server. Zhang et al. proposed a random projection ensemble approach DrugRPE to predict DTIs (Zhang et al., 2017), and several random projections build an ensemble REPTress system. In general, the method of fusing multiple predictors outperforms a single predictor. To facilitate comparison, the scores of accuracies (defined in Eq. (18)) obtained by these three predictors (He et al., 2010; Hu et al., 2016; Zhang et al., 2017) based on the benchmark datasets used in He et al. (He et al., 2010) via the 10-fold cross-validation test were listed in Table 2. Comprehensively, the comparative results showed that our model is more accurate than other existing methods.

GPCRs have proved to be one of the most important target families of modern drugs. Identifying the GPR-drug interaction is an important issue in bioinformatics, and a number of researchers have proposed effective predicted methods to identify GPCR-drupredictedions. Our method was also compared with the performance of different methods which predicting GPCR-drug interaction on the training dataset S G P C R − D r u g   over leave-one-out cross-validation, and validated in independent test dataset check390 (Xiao et al., 2013; Hu et al., 2016; Wang et al., 2020; Qiu et al., 2021). The results of the different methods tested on S G P C R − D r u g over leave-one-out cross-validation were shown in Table 3. The results of the other eight methods were reported in (Qiu et al., 2021). From Table 3, we can find that the MCC values of our method were 10% higher than others.

The generalization ability of machine learning models is usually evaluated through an independent test. The D92M is the GPCR-drug interaction dataset in (Wang et al., 2020), which is applied as a training dataset, and check390 as a validation dataset. The results of the validation test on check390 were listed in Table 4, which demonstrated that our method almost outperform the others across the five metrics, except for BOW-GBDT achieves the highest value of Sp (93.1%). Compared with other state-of-the-art methods, the ACC value of our method is 3.4% higher, the MCC value is 6% higher than the second one. All these results demonstrate the effectiveness of the proposed methods.