Image-based effective feature generation for protein structural class and ligand binding prediction

Structural class prediction dataset

We have used 40 percent ID filtered subset of PDB-style files for SCOPe domains version 2.03 (Fox, Brenner & Chandonia, 2013) as our dataset. It contains a total of 12,119 PDB files. Each PDB files contains SCOP(e) concise classification string (sccs) which respectively describes class, fold, superfamily, and family. In this paper, we are going to experiment only with the class of the protein. In the dataset, there are total seven protein structural classes. For benchmark analysis with CoMOGrad and Phog, the common PDB files were used as dataset. The common PDB files are total of 11,052. The details of the protein structural classes are given in Table 1. This dataset is widely used as a benchmark in the literature for protein structural similarity prediction (Karim et al., 2015).

Image generation

We have generated images of protein structures according to the methodology described in CoMOGrad and PHOG (Karim et al., 2015). Only α carbons of the amino acids in the protein structure are considered for image generation. From the three dimensional coordinates of the α carbon atoms, a distance matrix is generated by taking the Euclidean distance among all pairs. This distance matrix converts a 3D structure of a protein to a 2D matrix. Euclidean distance is applied as the distance measure because it exceeds the popularly applied costly alignment distance measure of α carbon distance matrices. Thus only half of the matrix contains redundant information due to symmetry. The matrix can be further regarded as an image.

Scaling of images

The dimension of protein images is based on the total number of α carbon they have. So, every individual protein images are different from the other in dimension. Therefore, the images were scaled to the same dimension. CoMOGrad and PHOG have used Bi-cubic interpolation and wavelet transform to scale all the protein images into 128 × 128 dimension (Karim et al., 2015). During the Bi-cubic interpolation step, most of the images were in 128 × 128 dimension so in the wavelet transform step they scaled all the images to that dimension. Thus, we have directly scaled the images to 128 × 128 dimension. We have used both real and scaled images to examine the differences in their predictive power. Sample rescaled images of protein structures are given in Fig. 2.

Feature extraction

We have generated five different feature groups. Our first four feature groups are different types of histograms and the fifth feature group is about the prognosis of the atoms. The histograms were taken from both scaled and unscaled images.

Local binary pattern histogram

Local binary pattern histogram was first proposed by Ojala, Pietikainen & Harwood (1994) and popularized by the work of Ojala, Pietikainen & Maenpaa (2002). LBP computes the local representation of the texture of an image as a texture descriptor. Comparing each pixel with its neighboring pixels the local representation is created. The image is transformed into a grayscale image. In a 3 × 3 neighborhood, the center pixel value is calculated by comparing with its eight neighboring pixels. Each comparison gives a result of either 0 if the center pixel value is greater than the comparing neighbor pixel or 1 for the latter. A clockwise direction starting from the top-left one provides a binary number. The binary number is converted to a decimal number and the value is placed in the center pixel. LBP codes or LBPs are the obtained binary numbers. An example of a basic LBP is given in Fig. 3. After calculating the value for each pixel of the image, a histogram is calculated. A 3 × 3 neighborhood has 2⁸ = 256 possible patterns, thus the values range from 0 to maximum 255 in each pixel of the image. The total number of bins of the histogram is thus 256. We would get 256 attributes from each image. We have used zero-padding technique to generate LBP.

Gabor filtered local binary pattern histogram

Gabor Filter is titled after Dennis Gabor. It is used for texture segmentation (Jain & Farrokhnia, 1991), optical character recognition (Jain & Bhattacharjee, 1992), edge detection (Mehrotra, Namuduri & Ranganathan, 1992), and so on. It is a linear filter that examines if there is any particular frequency content in the image in specific areas in a localized region throughout the point. In spatial and frequency domain, the Gabor filter has been determined to own optimal localization properties (Jain, Ratha & Lakshmanan, 1997). Gabor filter lets a particular band of frequencies while discarding others and that is why it is designated as a bandpass filter. It resembles the features of simple visual cortical cells (Dunn & Higgins, 1995). It provides the largest response where texture changes at edges and at points. The multiplication of a sinusoid and a Gaussian is called the Gabor filter. The filter has a real and an imaginary segment rendering orthogonal directions. The two segments can be set into a complex number or used individually. Various shapes, sizes and smoothness levels in an image can be detected by Gabor filter.

Equation 1 is the complex version of the Gabor function and Eqs. (2) and (3) are real and imaginary version respectively. The shape and size of the Gabor function is regulated by its five parameters. Here, λ controls the wavelength of this sinusoid, θ is the angle of the normal to the sinusoid, ϕ is the phase shift of the sinusoid, γ controls the aspect ratio, The spatial envelope or the standard deviation of the Gaussian is σ. For our experiments, we have used λ = 10, θ = 0, ϕ = 0, γ = 0.02 and σ = 5. After applying the Gabor filter, LBP techniques are applied to the image to get 256 attributes.

Atomic bond features

First of all, we’ve identified unique atoms amidst all the protein PDB files. From each protein PDB file, we’ve counted occurrences of each atom. Then we’ve taken the percentage as features of each atom among all the atoms that each protein has. Then we’ve taken the first 100 sequential atoms and used their atomic mass as the feature. Then we’ve counted the bond that each pair of atoms has in a particular protein using atomic distance based on a threshold value. Finally, we’ve taken the percentage as the feature of the bond of each unique pair of atoms among all the bonds that the protein has.

Separate row multiplication matrix with uniform LBP histogram

The image is split into 3 × 3 matrices. From each matrix, we get three rows with the dimension of 1 × 3. By multiplying each row with the same 3 × 3 matrix, we get three result matrix consisting of 1 × 3 dimension. Each cell is divided by 100. The results are then put in the 3 × 3 matrix in accordance with the row numbers. The color intensity of an image is between 0 to 255. So, if the value of any cell of the result matrix is greater than 255, then the value is replaced with 255. After applying this technique, the uniform LBP is applied. From Fig. 4, (a) presents a 3 × 3 section of matrix and the rows, (b) exhibits the result of multiplication, (c) shows the value after dividing by 100, (d) shows the replacement result of value greater than 255 and (e) shows a 3 × 3 matrix section after SRM-Matrix transformation.

Another variation of the LBP is called uniform pattern (Ojala, Pietikainen & Maenpaa, 2002). Some binary patterns occur more generally in texture images. If the binary pattern comprises at most two 0-1 or 1-0 transitions when the bit pattern is held circular then the pattern is called uniform. For instance, 01000000 has two transitions, 00000111 has two transitions which are uniform pattern on the other hand 01010100 has six transitions,11001001 has four transitions which are not uniform. A neighborhood with a dimension of 3 × 3 has 2⁸ = 256 possible patterns with 58 of them being uniform. For estimating the histogram, every uniform pattern gets a separate bin while a single bin is allotted for all non-uniform patterns. Therefore, from a uniform binary pattern, we get the histogram of the total bin size of 59.

Neighbor block subtraction matrix with uniform LBP histogram

Blocks are of the same dimension, 3 × 3. Two blocks of matrices are considered neighbors for this method if the center cells are neighboring. Because of this, the value of the last two columns of the first block and the first two columns of the second block are the same. The two blocks of matrices are subtracted and the result is set in the place of the first block. If any of the cells have any negative number, then 0 is placed instead of the negative value. The replacement of value is made because the histogram bin begins from zero. Uniform LBP is then used to compute the histogram.

Summary of all the feature groups used in this paper is given in Table 2.

Handling imbalance in data

From Table 1 it can be noted that the classes are imbalanced. To balance the classes, we have used SMOTE (Chawla et al., 2002). The percentage of SMOTE indicates that how many more instances would be generated. We have over-sampled our instances close to the highest number of instances among all the classes. If x denotes the highest number of instances among all the classes and y is denoted by a class which we will SMOTE, then the expression for the percentage calculation is ${\displaystyle \frac{x-y}{y}\ast 100}$. We have used five nearest neighbors to generate the over-sampled instances.

It is to be noted that SMOTE was used on only train dataset after separating the train and test dataset. So, there are no artificial instances in the test dataset.

Classifiers used

We have used five classifiers for the analysis of features applied to solve structural class prediction problem: K-Nearest Neighbor (KNN), Naive Bayesian Classifier, support vector machines (SVM), Adaptive Boosting (AdaBoost) and Random Forest. A concise description of the classifiers is given in this section.

K-nearest neighbor

K-nearest neighbor algorithm (KNN) (Aha, Kibler & Albert, 1991) is a similarity-based classification technique. It is a lazy classification technique. Distance metrics are used for each instance of the whole dataset for calculating the KNN. The labels of the nearest neighbors decide the label of the test instances. It works poorly for high dimensional data. Euclidean distance, Hamming distance, Manhattan distance, Minkowski distance, Tanimoto distance and Jaccard distance are used for similarity measures.

Naive Bayesian classifier

Naive Bayesian classifier (Maron, 1961) is based on probabilistic inference of samples observed where the decision variable and the features form a very naive structure of Bayesian Network. Naive Bayesian classifiers work best for image recognition and text mining.

Support vector machine

Support vector machine (Cortes & Vapnik, 1995) works by creating and separating hyperplane for a given dataset by sampling different classes which are separated by maximum width.

Adaptive boosting

Adaptive Boosting classifier (Freund & Schapire, 1997) is a meta-classifier which aims to make a strong classifier using a set of weak classifiers. The classifiers whose performance is marginally better than random classifiers are called weak classifiers.

Random forest

Random Forest (Ho, 1995) is an ensemble classifier. A decision tree is created in each iteration with features taken randomly. It samples selected features using bootstrap aggregating.

Ligand-binding dataset

We’ve used Computer Vision and Pattern Discovery for Bioimages Group @ BII as our dataset. In our dataset, 3,000 protein-ligand complexes were determined experimentally with 3D structures available. Each protein and its ligand are of one-to-one correspondence, that is, they can bind to each other and make Protein-Ligand complex. The dataset has 3,000 pairs of protein and ligand where the same name/ID of protein and ligand interacts/binds with each other.

We’ve used OpenCV (Bradski & Kaehler, 2008) library to create images from PDB files. For protein, we’ve considered the coordinates of only the alpha-carbons to generate the distance matrix to create an image. Because alpha-carbon can represent the structural information of protein quite well. But the given ligands were small in terms of atom number. So, while creating ligand images, we’ve considered all the atom’s coordinates for generating distance matrix.

Among the PDB files, 33 ligands have only one atom, which will create a 1 × 1 image having no significance for feature extraction. So, we had to compromise those 33 ligands as well as 33 corresponding proteins from training. Figure 5 shows a sample non-scaled and scaled image of a ligand.

Handling imbalance

The given dataset has only positive instances (the pairs of protein and ligand where they bind with each other). But there were no negative instances (the pairs of protein and ligand where they do not bind with each other). The missing negative instances have created our dataset highly imbalanced. To overcome this imbalance, we’ve generated negative instances in two different ways.

1. Random Negative Undersampling: We have 2,967 protein PDB and 2,967 ligand PDB where 8,803,089 pairs are possible. Among these, 2,967 pairs are given as positive instances and the rest 8,800,122 pairs are unknown/unseen instances. From the unseen pairs, we’ve taken 2,967 pairs randomly as negative instances to make our dataset balanced.

2. Clustering-Based Undersampling: Using the positive instances (2,967 pairs), we’ve created 10 clusters. Then we’ve searched for 2,967 unseen pairs randomly as negative instances where they belong to those 10 clusters. We’ve made sure that each cluster has the same number of positive and negative instances to make the dataset balanced (See Fig. 6).

Similarity-based classifier

We’ve developed a similarity-based clustering method to predict the binding class. Distance is used to measure similarity. Our methodology is given in Fig. 7 and the pseudo-code in Algorithm 1.

From the PDB dataset of proteins and ligands, firstly we have generated images and converted to 128 × 128 images for each protein and ligand. From these images, we have generated two different features.

1. CoMOGrad and PHOG: CoMOGrad stands for Co-occurrence Matrix of the Oriented Gradient of Distance Matrices and PHOG stands for Pyramid Histogram of Oriented Gradient (Karim et al., 2015). This methodology also uses the α carbon distance matrix of protein. The dimension of all distance matrix is converted to 128 × 128. In CoMOGrad, the gradient angle and magnitude is computed from the distance matrix and the values are quantized. Quantization is a compressing technique which compresses a range of values to a single quantum value. In this methodology, the values are quantized to 16 bins which produce a co-occurrence matrix which is 16 × 16 matrix. The matrix is converted into a vector of size 256. Quadtree from the distance matrix is created with the desired level in PHOG. Gradient Oriented Histogram of each node is calculated with the preferred number of bins and bin size. In gradient oriented histogram an image is divided into small sub-images called cells and histogram of edge orientations are accumulated within the cell. The combined histogram entries are used as the feature vector describing the object. Total features which are the multiplication of total nodes and number of bins are incorporated in the vector with the size of the total number of features. The vector is normalized by dividing it with the sum of its components.

2. Hybrid Local Binary Pattern: LBP (Ojala, Pietikainen & Harwood, 1994) is a procedure of LBP histogram. We have used all the five feature groups described in the last section for structural class prediction problem.

Distance can only be calculated between proteins or between ligands. We’ve used KNN and Clustering method to calculate these distances.

1. RELATEDLIGANDS ($\mathbb{N}\mathbb{P}$): For a given protein, find K-nearest proteins. The ligands those bind with the above nearest proteins are the Related Ligands for the given protein (See Fig. 8).

2. RELATEDPROTEINS ($\mathbb{N}\mathbb{L}$): For a given ligand, find K-nearest ligands. The proteins those bind with the above nearest ligands are the Related Proteins for the given ligand (See Fig. 9).

To find the distances between pairs of ligands and proteins are calculated using Euclidean and Manhattan distances. Threshold is the boundary between similarity and dissimilarity in terms of distance. If distance is less than the threshold, then prediction is positive similarity, else the prediction is negative similarity. The threshold of each category of distances is the average distance based on the number of nearest neighbors. Measuring distance has been done in two ways. One method is to get the cluster mean of k-RelatedLigands/RelatedProteins, then measure the distance from the GivenProtein/GivenLigand. Another method is to measure distances between the GivenProtein/GivenLigand and k-RelatedLigands/RelatedProteins, then take the average of those distances as the final distance.

For a given pair of Protein and Ligand, we want to predict if the will bind with each other or not. For measuring distance d_l, from the given protein, we searched for k-nearest proteins and found the k related ligands accordingly. Then we’ve calculated the distance using the above-mentioned methods. Then we’ve taken the vote for the binding class by all categories of distances based on their thresholds. Then finally, we’ve used majority voting mechanism to predict the binding class.

Hyperparameters

There are some hyperparameters of our proposed method.

1. Number of nearest neighbors: Our algorithm’s prediction accuracy is highly dependent on the number of nearest neighbors for finding both RELATEDLIGANDS($\mathbb{N}\mathbb{P}$) and RELATEDPROTEINS ($\mathbb{N}\mathbb{L}$). We’ve used five and three (best fit) nearest neighbors in this experiment.

2. Threshold: This is the threshold of distance for determining whether two proteins or two ligands are similar or not. For a higher value of threshold, there is a higher possibility for our algorithm to predict positive binding class for the majority of the protein-ligand pairs. Also, the lower the threshold is, the higher is the possibility of negative binding class prediction. We have taken the average of distances among three nearest neighbors as our threshold for each category of the distances. In terms of setting the threshold values, distances of each training sample from the whole training data was measured. While doing so, negative pairs of protein-ligand were considered when positive pair of protein-ligand was given for getting distances. Similarly, positive pairs of protein-ligand were considered when negative pair of protein-ligand was given. The average of these distances was taken as threshold values.