EnzyNet: enzyme classification using 3D convolutional neural networks on spatial representation

Dataset

The dataset used in this study has been retrieved from the RCSB PDB (http://www.rcsb.org), which contained 63,558 enzymes as of mid-March of 2017. It has been randomly split in training and testing set with the proportions 80/20%. Also, 20% of the training set has been put aside for validation and was later used for model selection. The details of each of those sets have been summed up in Table 1.

Representation space and occupancy grid

It has been established (Illergard, Ardell & Elofsson, 2009) that structure is far more conserved than sequence in nature. Since evolutionary relationships have been recognized as a confounding factor for EC, choosing an approach that focuses on spatial representation will be particularly robust in that respect. For this reason, we aim at building a model that seamlessly extracts relevant shape features from raw 3D structures. Accordingly, enzymes are represented as a binary volumetric shape with volume elements (voxels) fitted in a cube V of a fixed grid size l with respect to the three dimensions. Continuity between the voxels is achieved by nearest neighbor interpolation, such that for (i, j, k) ∈ [[0,l−1]]³ a voxel of vertices $$(i+\delta x,j+\delta y,k+\delta z)|(\delta x,\delta y,\delta z)\in {\{0,1\}}^3$$takes the value 1 if the backbone of the enzyme passes through the voxel, and 0 otherwise. Instead of extracting features from this 3D structure and analyzing the calculated multi-channel 2D feature maps (as performed in Zacharaki, 2017), we introduce this structure to the CNN, which optimizes its internal parameters (weights and biases of convolutional filters) using the backpropagation algorithm (Rumelhart, Hinton & Williams, 1986). The output of the convolutional filters comprises the feature maps that are selected by the method as best performing for the specific classification task.

To construct this shape representation, some preprocessing steps are necessary. First, protein structure is mapped to a grid of a predefined resolution. The selection of grid resolution determines the level of complexity/scale retained for the enzymatic structures. A full resolution is not preferred due to high data dimensionality and because fine local details are less relevant in characterizing enzymes’ chemical reactions. Thus, in order to avoid to get trapped into local minima, side chains are ignored and enzymes are represented exclusively through their “backbone” atoms that are carbon, nitrogen, and calcium.

Additionally, we note that enzymes do not possess any absolute spatial orientation. Unlike objects such as chairs or boats that appear usually with a specific orientation, proteins can have any orientation in 3D conformational space, thus the Cartesian coordinates defined by the model stored in PDB represent only a frozen in space and time snapshot of an overall highly dynamic structural diversity. The orientation is irrelevant to the properties of the protein. This observation underlines the need of either a rotation invariant representation or of a convention that makes output structures comparable one to another based on the definition of an intrinsic coordinate system. We define as origin of this intrinsic coordinate system the consensus barycenter of the protein as it is defined by taking into account only the four atoms of the backbone for each residue, and as axes the principal directions of each enzyme calculated by principal component analysis. Each structure is rotated around its center and the three principal directions of the enzyme aligned with the three axes of the Cartesian coordinate system. Instead of defining a common reference frame and aligning the objects before building the prediction model, other works (Brock et al., 2016; Hegde & Zadeh, 2016; Maturana & Scherer, 2015) applied rotations around relevant axes for data augmentation. We decided to spatially normalize the data instead of arbitrarily augmenting them (by random rotations) because the number of samples is already big enough and an adequate sampling of all possible orientations would lead to an extremely large dataset difficult to handle. However, similarly to these works the definition of orientation includes uncertainties about the direction (left-right, bottom-up), which we tackled also by data augmentation.

Another critical part of the process is to determine how the enzymes should be fit in their volumes, as those can be of all types of shapes and sizes. Should we scale enzymes separately to make them fit in their respective volumes? Or, on the contrary, should we scale all enzymes in a uniform manner? We choose to select the second option because of two reasons. First, doing otherwise would lead enzymes to be represented at different resolutions. Second, biological considerations invite us to make the convolutional network aware of the size difference between samples, as those may be an implicit feature regarding class determination. We already know that our source files provide the coordinates of the enzymes at a same scale. After all, proteins are comprised of various combinations of equally sized amino acids. This scaling issue is therefore equivalent to determining a maximum radius R_max so that the atom occupancy information contained in the sphere centered on the barycenter of the enzyme and of radius R_max fits into V. This situation is illustrated in Fig. 1.

R_max has to be large enough so that a sufficient number of enzymes fit the most of their volume inside V. Conversely, it also has to be small enough so that most enzymes are represented at a satisfactory resolution.

As a result, a homothetic transformation with center S and ratio λ defined by $$S\text{\hspace{0.17em}center\hspace{0.17em}of\hspace{0.17em}}V\text{\hspace{0.17em}and\hspace{0.17em}}\lambda =\lfloor \frac{l}{2}-1\rfloor \times \frac{1}{R_{\max }}$$ (1)is performed on all enzymes to scale them to the desired size.

When l is low, the grid is coarse enough so that the voxels of the structure have a contiguous shape. Conversely, big volumes tend to separate voxels, which engender “holes.” In that case, consecutive backbone atoms $(\overrightarrow{A_i},\overrightarrow{A_{i+1}})$ are interpolated by p regularly spaced new points computed by $$\frac{(p-k+1)\times \overrightarrow{A_i}+k\times \overrightarrow{A_{i+1}}}{p+1}$$ (2)where k varies from 1 to p. The latter is determined empirically beforehand on enzymes of the training set.

A last preprocessing step is to remove potential outliers from the volume. This is done by eliminating voxels that do not have any immediate neighbor.

The illustration of the output volume obtained for different grid sizes has been provided in Fig. 2.

Data augmentation

As noted earlier, differences in spatial orientation had to be considered so that volumes can be comparable one to another.

Keeping the orientation constant, we perform data augmentation by applying transformations that preserve the principal components along the three axes, i.e., flips and combination of flips. The number of possible transformations applicable to each protein is 2³−1.

During training, all samples go through the process described in Algorithm 1.

Architecture

We considered shallow architectures in order to develop a framework that can be trained with common computational means. A grid search of configurations has been conducted on the training set and led us to consider the two-layer architecture presented in Fig. 3: input volumes of size 32 × 32 × 32, containing the structural coordinates of the enzymes encoded in voxels, first go through a convolutional layer of 32 filters of size 9 × 9 × 9 with stride 2. Then, a second convolutional layer of 64 filters of size 5 × 5 × 5 with stride 1 is used, followed by a max-pooling layer of size 2 × 2 × 2 with stride 2. Finally, there are two fully-connected layers of 128 and six (the number of classes) hidden units respectively, concluded by a softmax layer that outputs class probabilities corresponding to the six first-level EC numbers.

Several components of the VoxNet architecture depicted by Maturana & Scherer (2015) have also been investigated in our network. Leaky ReLU with parameter α = 0.1 is used as activation layer after each convolutional layer. The L2 regularization technique of strength λ = 0.001 is applied on the network’s layers. Overfitting is also tackled using dropout throughout the network.

In summary, this model contains 804,614 distinct parameters (including biases), which is approximately 13% less than for the VoxNet architecture.

Additionally, the Adam optimizer introduced by Kingma & Ba (2015) has been chosen for its computational efficiency and low need of hyper-parameter tuning.

The categorical cross entropy is used as loss function, since the latter can be adapted to multiclass classification. Two approaches are considered for computing this loss. The first one assesses misclassification error irrespectively of individual class sizes, whereas the second one takes class imbalance into account and uses customized penalization weights for each class. In this second approach, the aim is to compensate the lack of training samples in under-represented classes by providing a greater penalization to the loss in the event of misclassification than for larger classes.

The cross entropy loss is given by $$ℒ=-{\displaystyle \sum _{x\in \text{set}}{\displaystyle \sum _{i=1}^6{w}_i}}·{\text{δ}}_{x,i}·\log ({\widehat{p}}_{x,i})$$ (3)where ${\widehat{p}}_{x,i}$ is the predicted probability of enzyme x belonging to class ECi, δ_x,i the quantity equal to 1 only if enzyme x belongs to class ECi, and w_i the weight associated to class ECi. For the first approach, all weights w_i are taken equal to 1. For the weighted approach, we chose weights according to the formula $$w_i=\frac{\underset{j\in [\text{}[1;6]\text{}]}{\max }\#(\text{EC}j\text{\hspace{0.17em}training enzymes})}{\#(\text{EC}i\text{\hspace{0.17em}training enzymes})}$$ (4)which increases the contribution of the under-represented classes by an amount inversely proportional to their size and in respect to the largest class.

Metrics

Among the various multi-class metrics that have been studied by Sokolova & Lapalme (2009), we selected the most representative ones to assess the model’s performance. The metrics are based on the confusion matrix whose elements C(i, j) with i, j ∈ [[1, 6]] indicate the number of enzymes that belong to class ECi and are predicted as belonging to class ECj. They include:

• Accuracy which captures the average per-class effectiveness of the classifier: $$\text{Accuracy}=\frac{{\displaystyle \sum _{i=1}^{6}C}(i,i)}{{\displaystyle \sum _{i,j=1}^{6}C}(i,j)}$$

• Precision, recall and F1 score which are calculated per class ECi: $${\text{Precision}}_{\text{EC}i}=\frac{C(i,i)}{{\displaystyle \sum _{j=1}^{6}C}(j,i)}{\text{Recall}}_{\text{EC}i}=\frac{C(i,i)}{{\displaystyle \sum _{j=1}^{6}C}(i,j)}$$ $$F{1}_{\text{EC}i}=2\times \frac{{\text{Precision}}_{\text{EC}i}\times {\text{Recall}}_{\text{EC}i}}{{\text{Precision}}_{\text{EC}i}+{\text{Recall}}_{\text{EC}i}}$$

For each class, precision gives us an idea of the proportion of correctly classified enzymes among enzymes that have been classified in that class, while recall highlights the proportion of correctly classified enzymes among enzymes that actually belong to that class.

• Macro precision, recall and F1 score which express average performance over the six enzyme classes: $${\text{Precision}}_M=\frac{1}{6}{\displaystyle \sum _{i=1}^6{\text{Precision}}_{\text{EC}i}}{\text{Recall}}_M=\frac{1}{6}{\displaystyle \sum _{i=1}^6{\text{Recall}}_{\text{EC}i}}$$ $$F{1}_M=2\times \frac{{\text{Precision}}_M\times {\text{Recall}}_M}{{\text{Precision}}_M+{\text{Recall}}_M}$$

Final decision rule

At testing time, several approaches are considered for determining an enzyme’s final class.

• The “None” strategy consists of making a prediction based on the model of the 3D shape without any transformation.

• In the “flips” approach, the 2³−1 possible combinations of flipped volumes are generated and introduced to the classifier. The final class prediction $\widehat{i}$ is either:

  • –

    the class of maximum total probability (probability-based decision), determined by $$\widehat{i}=\arg \underset{i\in [\text{}[1,6]\text{}]}{\max }{\displaystyle \sum _{j=1}^n{\widehat{p}}_{i,j}}$$with n the total number of flips used during inference and ${\widehat{p}}_{i,j}$ the probability to belong to class i during flip j

  • –

    or the class selected by majority voting (class-based decision), computed through the formula $$\widehat{i}=\arg \underset{{}_{i\in [\text{}[1,6]\text{}]}}{\max }{\displaystyle \sum _{j=1}^n{\widehat{\delta }}_{i,j}}$$where ${\widehat{\delta }}_{i,j}$ equals 1 if at flip j, the enzyme is predicted to belong to class i, 0 otherwise.

• The “weighted flips” (W. flips) strategy is also based on fusion of decisions produced for each flipped volume, but this one weights each decision by a different coefficient, such as $\frac{1}{{\text{δ}}_x+{\text{δ}}_y+{\text{δ}}_z+1}$, which highlights transformations with the least number of flips. Just like the previous approach, the final result can be determined from either a probability- or class-based viewpoint.

Hyperparameter selection

Radius R_max and interpolation parameter p

A crucial point in the presented approach is to make a cogent selection of R_max. As previously discussed, this parameter controls the trade-off between the level of information (l) retained in each volume and the resolution with which they are conveyed (R_max). Figure 4 shows the analysis of the dataset from these two perspectives. The graph on the left helps us assess the minimum radius for which a decent amount of enzymes will be totally included in the volume, whereas the graph on the right highlights the quantity of information retained by each radius.

In details, from Fig. 4A we can see that the majority of enzymes can fit in a sphere with radius between 25 and 75 Å, with a peak around 35 Å. In Fig. 4B, each point (x, y) on a curve of radius R shows the percentage y of training enzymes that have at least x points included within radius R around their respective enzyme barycenter. Based on both graphs, we select a value of R_max = 40, that is big enough to capture more than half of enzymes in V, but small enough so that the smallest enzymes have radii of at least half of R_max.

Empirical observations show that a grid of l = 32 does not require any atom interpolation. Therefore, p is set to zero for our computations. For denser grids, e.g., with grid size l = 64 or 96, appropriate p values were p = 5 and p = 9, respectively.

Computation time

Our architecture was implemented on Python 3 using Keras (Chollet, 2015) on top of a GPU-enabled version of TensorFlow (Abadi et al., 2015). Enzyme information has been extracted using the open-source module BioPython, a fast and easy-to-use tool presented by Cock et al. (2009). The complete training of our model on 50,846 samples took about 5 h on an Intel i7 6700K machine with 32 GB of RAM and a GTX 1080 graphics card. During the training phase, the code computes batches in parallel and stores them in a queue in real-time, so that the main computational burden comes from neural network calculations and not the representation process. The average prediction time of the function of a single enzyme was about 6 ms without flips, or 50 ms with flips.

Further improvements of the method

In parallel, we investigated possible improvements in the architecture. We introduced batch normalization (Ioffe & Szegedy, 2015) at several positions of the network and repeated the same experiments as before. Those were placed after every activation layer, as suggested by Mishkin, Sergievskiy & Matas (2017). Also, Leaky ReLU activation layers were replaced by PReLU (He et al., 2015) layers which enable the network to learn adaptively the Leaky ReLU’s best α parameter. Although these changes led the network to converge faster to a stable optimum, the final performance increase was of an order of magnitude of only one percent and came at a higher computational cost.

In the following, we considered the transition from a binary representation capturing only shape information to “gray-level” representation capturing also information content. Several biological indicators such as the hydropathy index (Kyte & Doolittle, 1982), or isoelectric points (Wade, 2002) can be used to better describe local properties of amino acids that are the building blocks of the protein structure. From the perspective of computational analysis, these attributes can be incorporated into the representation model by attaching to the shape also appearance information. Figure 6 illustrates this idea showing the volumetric representation of two different attributes (hydropathy on the left and charge on the right). The different attributes, if up to three, can also be merged into a single structure visualized in color (RGB) scale. Different approaches to handle multi-channel volumetric data (images) by CNNs have recently been presented in works such as those of Zhang et al. (2017), Cao et al. (2016), Kamnitsas et al. (2017) and Nie et al. (2016). The literature however on the use of deep networks for 3D shapes with multi-channel appearance is limited. As preliminary analysis, we applied EnzyNet without modification on the architecture, but with re-training for the adjustments of the weights of the convolutional kernels. We introduced to the network either a single attribute (as a binary 3D image or gray-level 3D image), such as shape, hydropathy and isoelectric points, or a combination of these different attributes. Two other ways of information combination were examined: either each attribute was introduced as a different channel in the CNN, or the outputs of the single channel networks were combined through a fusion rule. However, all methods showed an increase of the order of magnitude of one percent. Further investigation will be required to appropriately harness this added information.

Furthermore, in the current study a relatively small grid size was selected due to computational limitations. However, as can be seen in Fig. 2 the details of the spatial structure of proteins can be better differentiated for higher values of l. The next step of this work would be to adapt the current architecture accordingly into a similar network that processes higher-resolution volumes. This could enable the network to capture more subtle features and potentially boost the performance of the classifier. On the other hand, the model performance using a finer representation might be more sensitive to possible inaccuracies in protein structure prediction. Initial tests showed that the appearance of artifacts in the location of atoms at l = 64 or 96 had a negligible impact on accuracy.

Wei et al. (2015) present a method that enables 3D CNN to deal with multi-label classification problems. This approach is interesting for us as it allows the extension of our method to the classification of multi-label enzymes. The obtained performance could subsequently be compared to previous work Amidi et al. (2017).

Other approaches based on the 3D representation of enzymes are also possible. Brock et al. (2016) performed very well on the ModelNet dataset using generative and discriminative modeling. Their voxel-based autoencoder is helpful for assessing the key features that are correctly learned from the 3D shapes. We could elaborate this information in order to identify significant features in a pre-training phase aiming to obtain better prediction performance.