Prediction of bacterial E3 ubiquitin ligase effectors using reduced amino acid peptide fingerprinting

Dataset

We identified a set of 164 confirmed bacterial or viral E3 ubiquitin ligase effectors from the UniProt database by searching for ‘E3 ligase’ in manually annotated bacterial and viral sequences and manually checking the results for accuracy (Bairoch et al., 2005). Negative examples were 235 other bacterial effectors identified from literature (Lee, Mazmanian & Schneewind, 2001; Stebbins & Galan, 2001; Bairoch et al., 2005; Burstein et al., 2009; Quezada et al., 2009; Samudrala, Heffron & McDermott, 2009; Spallek, Robatzek & Gohre, 2009; Buchko et al., 2010; Collins & Brown, 2010; Hicks & Galan, 2010; Price & Kwaik, 2010; Komander, Clague & Urbe, 2009; Wu et al., 2010; Dean, 2011; Lin et al., 2011; McDermott et al., 2011; Anderson & Frank, 2012; Deslandes & Rivas, 2012; Xin et al., 2012). Though somewhat limited in size, this set of E3 ubiquitin ligase virulence effectors represents the state of current knowledge, and we have had success with similar approaches applied to smaller datasets previously (Samudrala, Heffron & McDermott, 2009; McDermott et al., 2015). We include details on the dataset as Supplemental Data.

To provide predictions for relevant bacterial pathogens we downloaded a set of 171 genomes that are listed as human pathogens and are representative reference genomes from PATRIC (Wattam et al., 2017). This set comprises 480,562 protein sequences excluding all of the proteins used in the training set above. A list of the genomes included in this study is provided in Supplemental Data.

Features

Every protein sequence used for either learning or prediction is encoded by counting occurrences of peptides of varying length in the sequence in a manner similar to the previously described string kernel (Leslie, Eskin & Noble, 2002). The possible number of peptides greater than 4 amino acids long is very large (20⁴ = 160,000 peptides). We wanted to extend this approach to identify sequence patterns based on groupings of amino acids based on physiochemical or other properties. We therefore also encoded sequences to reduce the sequence space using one of several encodings (Table 1). Features were then generated for a range of different peptide lengths (3 to 20) and peptides that were observed in fewer than 10 examples were removed from consideration.

Features for each protein are generated by considering all peptides of length k in a sequence, encoding these (optionally) using the chosen encoding scheme, then counting the occurrences of the encoded peptide.

Data partitioning

To remove bias created by having multiple examples with very similar features (i.e., closely related effectors from different organisms) we first partitioned the examples to identify clusters of related effectors. In order to achieve this partitioning, we clustered the sequences based on NCBI BLASTP similarity results. Parameters of BLASTP were set to their default values. Using a lower E value threshold (for example, E = 0) groups sequences more tightly and thus results in clusters that are likely to be more similar to another cluster and thus represent a generous division of families for the classification task using our cross-validation approach (see below). Conversely, higher E value thresholds (for example, E = 0.01) yield broader, more general clusters that are less likely to be similar to any other clusters, and thus represent a conservative division of families for our classification task. We used a more conservative threshold (E 0.01) to group the set of 407 proteins into 176 clusters of loosely related protein sequences. We examine the effect of varying the BLAST E-value threshold on the size of the generated protein families (Fig. S1). We next compared this approach to a more sophisticated approach for determining orthologous groups of proteins OrthoMCL, and found only minor differences with our approach resulting in joining three of our previous clusters (each containing a single protein each) in to one. We then used this final set of 174 clusters for our analyses.

Cross validation

Cross validation (CV) is widely used to test the performance of a classification scheme on a given dataset. The entire dataset is partitioned into several non-overlapping folds. These folds are used as test sets. The corresponding training set for a particular fold consists of the remainder of the dataset. Each iteration of cross validation involves using a training set to generate a model and testing that model on the corresponding test set. This process is repeated until every fold has been tested.

The experimental setup of our study uses a variant of CV called Family-Wise Cross Validation (FWCV) to judge the performance of our classifier. FWCV places all the samples belonging to a particular cluster (see above) in a single test set, while the classifier is trained using the remaining data. This prevents model overfitting by reducing the trivial similarities between testing and training sets (i.e., those similarities based on traditional sequence similarity).

We use the following example to explain this process. In Table 2 protein sequences have been assigned a sequence family based on sequence similarity using traditional methods (like OrthoMCL). A FWCV run would select, for example, sequence families 1 and 4 to train on for a single fold. This would mean that sequences A, B, G, and H would be in the training set and sequences C, D, E, and F would be set aside for testing. A model would be trained on A, B, G, and H then applied to C, D, E, and F to assess performance. A good performance (as assessed by AUC) would mean that the information in sequences A and B from sequence family 1 could be used to predict the class of sequences C and D in sequence family 2. The process is then repeated for a number of folds (not useful for this limited example, but very useful in a real dataset), and the performance of the individual folds averaged to get an overall assessment of model performance.

Classification

The Support Vector Machine (SVM) determines the optimally separating hyperplane between two sets of points in high-dimensional feature space each belonging to a different class (Noble, 2006). We utilized the radial kernel from the e1071 R library (version 1.6-8; Meyer et al., 2017) in our implementation.

The area under the curve (AUC) and receiver-operator characteristic curve (ROC) calculation was performed using the R library pROC (version 1.10.0; Robin et al., 2011).

Feature selection

Feature selection was accomplished using two complementary methods. The first is the standard SVM Recursive Feature Extraction (SVM-RFE; algorithm as described in (Guyon et al., 2002)). We can obtain an ordering of the features using the absolute value of the entries of the SVM weight vector w. Each recursive feature elimination iteration involves eliminating the set of features that have the smallest absolute weight w_i until k features remain. The smaller and smaller sets of features obtained at each step represent predictive models that consider all the features together, but features selected vary due to the cross-validation approach we are using. In the second approach we established a simple metric to score individual features based on their representation in the positive example families versus negative example families. This was accomplished by calculating the score (S) for each feature (f) as the percentage of examples from each sequence cluster (M) that contained f, then calculating difference between the mean percentages for the set of positive families and negative families. Thus a positive score for a feature means that it is disproportionately represented in the positive examples while accounting for differences in the size of sequence families. Sets of individual predictive features were then used to train a minimal model for prediction as described in results.

Implementation details and availability

Feature generation from sequences is performed using a standalone Python (version 3.6.3; (Foundation, 2018) script that uses the BioPython library (version 1.70; (Cock et al., 2009)). Training and validation of models was performed in R (version 3.3.3 (R Core Team, 2017). The SVM-RFE algorithm used by SIEVE-Ub was implemented in R as described by GIST-RFE (Guyon et al., 2002).

Code for the algorithm and datasets used to produce the results described in this paper are available at https://doi.org/10.6084/m9.figshare.7766984.v1 and the current version of SIEVEUb is available at https://github.com/biodataganache/SIEVE-Ub.

Dissimilar ubiquitin ligases can be detected using reduced amino acid (RED) peptides

To provide feature sets that were specific enough to capture relationships between functionally similar proteins, yet general enough to identify regions of similarity between divergent sequences we adapted the kmer approach to represent protein sequences by a series of all peptides of length k from that sequence. Our novel extension translates each amino acid in the sequence to a smaller number of groups based on physicochemical properties or other arbitrary grouping methods- a reduced amino acid (RED) alphabet. Initially we chose three reduction mappings based on previously reported approaches: hydrophobicity (RED0), standard physiochemical properties (RED1), and solvent accessibility (RED2) (Arnold et al., 2009; Bacardit et al., 2009). The groups are listed in Table 1.

The set of positive and negative examples for E3 ubiquitin ligases was encoded using each of the REDs and the native sequence, and peptide kmers of various lengths were counted for each. Peptides present in fewer than 10 examples were excluded from further consideration. Each dataset was then split into independent training and testing sets on a sequence cluster-wise basis (that is, clusters of similar sequences as determined by OrthoMCL were kept together in the training or testing set), based on a conservative cluster grouping (see methods). Cluster-wise splits and associated training and testing were performed 100 times for each model and the score (SVM discriminant) for each example averaged. Average scores were used to determine ROC AUC for each model and results are presented in Table 3 and Fig. S1. We were concerned that the cross-validation procedure might be affected by the large size of some of the sequence families (an NEL cluster has close to 100 members) or by how the families were chosen during cross-validation. We examined the first problem by reducing the size of the large families to a maximum value and rerunning the analysis, finding that this did not affect the results. Similarly, choosing different random seeds prior to cross-validation did not significantly affect results (see Supplemental Data).

Surprisingly, the models using RED0, a simple division of amino acids into hydrophobic and hydrophilic residues, performed the best for nearly all peptide lengths with a maximum AUC of about 0.85. The maximum AUC observed occurs with RED0 and a peptide length of 14 (RED0-K14) and so we focused on characterization of this model for the remainder of the article. Our results indicate that a simple encoding of amino acids can be used to classify effectors with E3 ubiquitin ligase function from other effectors, and from other non-effector proteins in general (see Prediction of novel E3 ubiquitin ligase mimics, below), with good confidence.

We hypothesized that the performance of the RED0 is based on accurately representing the pattern of hydrophobic and hydrophilic residues in kmers. To examine this hypothesis we applied a family-wise cross-validation approach using ten alphabets where residues had been randomly assigned to either the hydrophobic or hydrophilic groups preserving the overall balance of hydrophobic to hydrophilic residues in the resulting random alphabet (6:14; see Table 1). We compared the performance of these random binary REDs at a kmer size of 14 with the true hydrophobic/hydrophilic RED0-K14 also run ten times to show the variability in partitioning of training and testing sets inherent in our approach and show the results in Fig. 1. In all cases the true RED0 outperforms the randomized REDs supporting our hypothesis though we note that there is a wide range of performances given with random binary REDs. We believe this is due to some random assortments containing reasonable divisions of residues between hydrophobic and hydrophilic residues because of the very simple nature of this division.

SIEVE-Ub identifies biologically functional peptides

To identify a minimal set of features that are important for classification of E3 ubiquitin ligases from other effectors we used recursive feature elimination, a standard machine learning approach (Samudrala, Heffron & McDermott, 2009). Briefly, a model is trained on all features using our family-wise cross-validation approach, then weights for each feature are used to discard 50% of the features with the lowest impact on model performance. The remaining features are then used in another model training round in which this process is repeated until all the features have been eliminated. Using this approach we found that performance of the model dropped off when the number of features was still quite large, >2,000, and so does not identify a minimal set of features important for discrimination of examples. The training performance results from the RFE on the RED0-K14 model are shown in Fig. S2.

Since RFE failed to identify a minimal set of predictive features we developed a simple scoring metric to evaluate each feature independently to identify those features with disproportionate representation in the positive example set, while accounting for differences in sequence family sizes (see Methods). We applied this score across various kmer lengths and REDs and show the results in Fig. 2. Similar to the results we obtained in performance of models incorporating all features (Fig. S2), this approach shows that RED0 results in the longest kmers that are specific to positive examples with lengths of 12–14. This also shows that the naturally occurring amino acid sequence does not produce kmers specific to the positive examples, highlighting the strength in using our RED approach.

We trained models using our family-wise cross-validation approach with the top most predictive kmer features from RED0 with a kmer length of 14 for consistency with our previous results, and found that the most predictive model performed quite well (AUC 0.87) with just ten features. The features from the minimal model are provided as Supplemental Data along with their locations in each of the positive and negative examples in our analysis set.

Though the E3 ligase examples used as our positive examples are diverse in terms of sequence many do fall into the families of E3 ligases described in the Introduction; HECT/U-box, RING, and NEL. We hypothesized that the most predictive kmers identified in our analysis would map to the known E3 ligase domains in these proteins. We show in Fig. 3 our scoring metric across four examples from different E3 ligase families. Significantly predictive peaks corresponding to highly predictive kmers can be found in each of the known domains from the HECT-like, NEL, AvrPtoB, and recently identified RavN effectors demonstrating that our approach is able to establish connections across disparate sequence families. We note that, except in the case of AvrPtoB, the most predictive kmers map to regions outside the known E3 ligase domains, indicating the presence of other signals that might be important in prediction.

We also applied our model to the recently discovered E3 ubiquitin ligase effector RavN from Legionella, which was not included in our training set (Lin et al., 2018). Our method predicted RavN to be an E3 ubiquitin ligase with a probability of 85%, despite having no detectable sequence similarity with other known E3 ligases. Though this is a limited validation it demonstrates the power of our approach at identification of novel E3 ligase effectors that have vastly divergent sequence.

Prediction of novel E3 ubiquitin ligase mimics

To predict novel E3 ubiquitin ligase mimics in a larger set of sequences we applied the model described above (kmer 14 in RAA0, top 2000 most important features) to a set of over 400,000 proteins from representative human pathogens obtained from the PATRIC database (Wattam et al., 2017). We found only 67 proteins with positive SIEVEUb scores greater than 0.5, indicating that potential ubiquitin ligases are not common as assessed by our approach. We note that PATRIC annotations include E3 ubiquitin ligase functions, and none are detected in the subset we’ve focused on. This indicates that the predictions we’ve made are truly novel. We further filtered this list to include predictions that occurred in bacteria containing type III, IV, or VI secretion systems. This yields a list (Table 4) of predicted E3 ubiquitin ligases that are in organisms capable of delivering effector proteins in to the eukaryotic host cell cytoplasm, though we note that such effectors could be secreted via other mechanisms. Several of these predictions are annotated as enzymes, which could be false positive predictions. However, many virulence effectors are known to be multifunctional and annotation of functions is prone to error.