Squidly harnesses enzyme functional hierarchy and contrastive learning to efficiently predict catalytic residues from sequence
Figures
Overview of the model.
(A) Depicting a sequence being translated into per-token embeddings via ESM2, each per-token embedding is then translated into a smaller representation space by the MLP trained with contrastive loss (CL). Each residue from the CL model is then predicted as catalytic or not catalytic via a bidirectional LSTM network. The CL and LSTM models are ensembled (5 models) to improve performance and capture variance between the models. (B) Representation of the pair-mining schemes depict the negative and positive classes for each scheme. Scheme 1 is the most permissive and is amino acid and EC unaware. Schemes 2 and 3 are more restrictive, requiring both the amino acids and the EC classes to be shared or different in the pair-mining scheme. Scheme 3 included additional negative pairs to separate EC numbers in the latent space.
Performance of Squidly pair-mining schemes on the Uni3175 dataset.
(A) We compare the Squidly schemes without the ensemble with AEGAN, the tool by Shen et al., the authors of this dataset (Shen et al., 2023). Reaction informed schemes (Schemes 2 and 3) achieved the highest F1 scores, slightly surpassing AEGAN when using the 15B ESM2 model. Results are reported as the mean across five independent runs (n=5); error bars indicate standard deviation. The 3B ESM2 model remains competitive at a much lower computation cost. In comparison to the rationally informed pair schemes, the Scheme 1 (random pairing, S1) and raw embedding LSTM models performed poorly. (B) Sequence similarity based on sequence identity to the closest sequence in the training set between each family. (C) F1, recall, and precision for the ensemble of Squidly (3B, scheme 3) and BLAST with different sequence identity cut-offs. The cut-offs indicate the point at which to transition from a Squidly to a BLAST prediction. At 100% only Squidly is used for predictions, and at 0%, BLAST is used, unless there are no similar sequences. For sequences with less than 30% sequence identity (dashed black line), Squidly predictions are used. The dashed horizontal line represents the total score for BLAST on the specific dataset, while the dots represent the combined BLAST and Squidly ensemble approach. Squidly reports higher scores for recall, while BLAST is more precise. The performance is best when used as an ensemble.
The default ensemble threshold for prediction and variance cutoffs of 0.6 and 0.225 (outlined in red) were selected to balance the precision and recall of the Squidly ensemble model on the Uni3175 benchmark under varying mean prediction and variance thresholds.
Shown are F1 scores, precision, and recall for the 3B and 15B models. (A) F1 of the 3B model, peaking at 0.70. (B) F1 of the 15B model, peaking at 0.69. (C) Precision of the 3B model. (D) Precision of the 15B model. (E) Recall of the 3B model. (F) Recall of the 15B model.
Squidly performance on the CataloDB catalytic-residue prediction benchmark.
(A) BLAST sequence similarity to the training set for low identity CataloDB test set sequences. (B) Structural similarity for the final test set, filtered on both structure and sequence. (C) F1, recall, and precision for Squidly (scheme 3), SCREEN and BLAST on CataloDB. (D) The number of CataloDB catalytic residues predicted (true positives) by Squidly (scheme 3) for the six available EC classes relative to the decreasing total number of true catalytic residues. (E) The precision, F1, and recall score of Squidly predictions for the six available EC classes. (F) The distribution of the identity of the catalytic residues is not even and represents the uneven distribution of the catalytic mechanism within the benchmark dataset. (G) Almost 50% of the data points are in the well characterised class of EC 3, which are hydrolases. EC classes 6 and 7 are poorly represented: ligases (EC 6) have only three representatives, while translocases (EC 7) are absent entirely. (H) Performance of the Squidly 15B ensemble model on the CataloDB benchmark under varying prediction and variance thresholds. Note, default 0.5 prediction cutoffs were used to evaluate Squidly’s performance in (C).
Performance of the Squidly ensemble model on the CataloDB benchmark under varying mean prediction and variance thresholds.
Shown are F1 scores, precision, and recall for the 3B and 15B models. (A) F1 of the 3B model, peaking at 0.70. (B) F1 of the 15B model, peaking at 0.69. (C) Precision of the 3B model shows that as the cutoff increases. (D) Precision of the 15B model. (E) Recall of the 3B model. (F) Recall of the 15B model. The default ensemble threshold for prediction and variance cutoffs of 0.6 and 0.225 (outlined in red) were selected using the Uni3175 benchmark.
The maximum F1 score achievable with the Squidly ensemble model for each EC number in the CataloDB benchmark under varying mean prediction and variance thresholds.
(A) F1 score as a function of the mean prediction threshold, with the variance threshold fixed at 0.225. (B) F1 score as a function of the variance threshold, with the mean prediction threshold fixed at 0.6. The default threshold of 0.6 and 0.225 determined using the Uni3175 benchmark is relatively stable for each EC number. However, the limited data for EC 5 and 6 (N=8, N=3, respectively) is such that users are recommended to further test the decision thresholds for their intended application.
Amino acid distribution of catalytic residues in dataset 1.
Dataset 1 contains 3602 catalytic residues, with 15 possible amino acids. The database is skewed towards certain amino acids which are commonly involved in catalysis.
Dataset 1 EC distribution.
Dataset 1 is made up of 2030 sequences from Swissprot. The sequences have experimental validation for the catalytic residue annotations. Validation sets are derived exclusively from this distribution. A clear bias exists for hydrolases (EC 3 .X).
Amino acid distribution of catalytic residues in dataset 2.
Catalytic residue distribution. Dataset 2 contains 9781 catalytic residues, with 18 possible amino acids. Amino acids M, V, and F are additions not seen in dataset 1, although they exist in very few numbers.
Dataset 2 EC distribution.
Dataset 2 is made up of 5921 sequences. The sequences are taken from all the available proteins with catalytic site annotations in Swissprot. A very similar distribution is seen between datasets 3 and 2, with a notable increase in EC 2 again.
Amino acid distribution of catalytic residues in dataset 3.
Dataset 3 contains 44,252 catalytic residues, with 19 possible catalytic amino acids. Isoleucine is an additional amino acid not seen in datasets 1 or 2. The distribution seen here is very similar to that of datasets 1 and 2.
Dataset 3 EC distribution.
Dataset 3 is made up of 27,401 sequences. The sequences are taken from all the available proteins with catalytic site annotations in Swissprot. A very similar distribution is seen between datasets 3 and 2, with a notable increase in EC 2 again.
Principal component analysis of per-residue embeddings from EC numbers 3.1 and 2.7.
In the upper row, red colours indicate a catalytic residue, and blue, non-catalytic. In the lower row, residues are coloured by amino acid class which represent their structural and chemical properties. The x and y axes represent the first and second principal components in each subplot. A clear separation of catalytic and non-catalytic residues is seen in the principal components which explain the most variance in our data.
Catalytic-residue prediction by Squidly, EasIFA, and AEGAN on the filtered EasIFA benchmark.
(A) Recall, precision, and false positive rate (FPR) for AEGAN, EasIFA, and Squidly on the filtered EasIFA test subset. EasIFA shows balanced performance, Squidly is competitive given its sequence-only inputs, and AEGAN displays unusually low precision that may reflect a benchmarking artefact. (B) Most of the test set are sequences from EC classes 2 and 3, with no representation for classes 4, 6, and 7. (C) The differences in false positives and true positives predicted by Squidly’s best model and EasIFA. Although the models agree with 100% of true positive predictions made by both tools, Squidly has a higher tendency for false positives, with many of the false positives being true binding sites.
Tables
F1 scores of tools on six common benchmark datasets.
| Tool | EF fold | EF superfamily | EF family | HA superfamily | NN | PC |
|---|---|---|---|---|---|---|
| PREvaIL (Song et al., 2018) | 0.2579 | 0.2575 | 0.2585 | 0.2483 | 0.2537 | 0.2582 |
| CRpred (Zhang et al., 2008) | 0.2613 | 0.2643 | 0.2651 | 0.2628 | 0.2638 | 0.2630 |
| AEGAN* (Shen et al., 2023) | 0.8369 | 0.8687 | 0.8552 | 0.8259 | 0.8177 | 0.8736 |
| SCREEN* (Pan et al., 2024) | 0.645 | 0.615 | - | 0.72 | 0.738 | 0.741 |
| Squidly (ESM2 3B, S3) | 0.8631 | 0.8775 | 0.8545 | 0.8171 | 0.8354 | 0.8527 |
| Squidly (ESM2 15B, S3) | 0.8765 | 0.8936 | 0.8867 | 0.8513 | 0.8586 | 0.8454 |
| BLAST (ultra-sensitive) | 0.9437 | 0.9216 | 0.8967 | 0.9053 | 0.9169 | 0.9119 |
| BLAST +Squidly (ESM2 3B, S3) | 0.9509 | 0.9333 | 0.9126 | 0.9111 | 0.9259 | 0.9119 |
| BLAST +Squidly (ESM2 15B, S3) | 0.9476 | 0.9274 | 0.9034 | 0.9179 | 0.9230 | 0.9119 |
-
Benchmark results for PREvaIL, CRpred, and AEGAN were sourced from Shen et al., 2023. Results from Screen were sourced from Pan et al., 2024. BLAST results were performed using the ultra-sensitive flag and aligned to the closest recorded enzyme with an active site in the training dataset (Uni14230). Squidly results shown for models trained on the Uni14230 dataset (same as AEGAN) using pair mining scheme 3 (S3).
-
*
Indicates structure-based tools.
Resource usage measurements for Squidly models compared to structure prediction tool Chai-1 when predicting 100 sequences with an average length of 443 residues.
| Model | Wall clock (s) | VRAM (GB) | RAM (GB) |
|---|---|---|---|
| Squidly 3B | 43.9 | 26.5 | 17.1 |
| Squidly 15B | 185.6 | 74.8 | 91.5 |
| Chai-1 | 13263.0 | 20.8 | 19.3 |
-
All models benchmarked using the Nvidia H100 GPU with 80 Gb of VRAM, AMD EPYC 7643 48-Core Processor on a HPC cluster running Rocky Linux 8.10.
-
Squidly models are ensembled (n=5).
Benchmark datasets in prior work.
| Dataset | Year of publication | Size (before sub-setting) | SCREEN subset | SCREEN >90% identity | AEGAN subset | AEGAN >90% identity |
|---|---|---|---|---|---|---|
| EF family (Youn et al., 2007) | 2007 | 293 | NA | NA | 196 | 3.6% |
| EF superfamily (Youn et al., 2007) | 2007 | 189 | 49 | 14.3% | 93 | 2.2% |
| EF fold (Youn et al., 2007) | 2007 | 236 | 60 | 10.0% | 123 | 3.3% |
| HA superfamily (Chea and Livesay, 2007) | 2007 | 282 | 191 | 77.0% | 152 | 2.6% |
| NN (Bartlett et al., 2002) | 2002 | 163 | 95 | 83.1% | 108 | 2.7% |
| PC (Petrova and Wu, 2006) | 2006 | 81 | 48 | 85.4% | 55 | 1.8% |
| SCREEN (Pan et al., 2024) training | 2024 | 1055 | 1055 | NA | NA | NA |
| AEGAN (Shen et al., 2023) training (Uni14230) | 2023 | 8435 | NA | NA | 8435 | NA |
| AEGAN test (Uni3175) | 2023 | 1955 | NA | NA | 1955 | 0.00% |
Brief description of methods compared to Squidly in common benchmarks.
| Study | Input | Method |
|---|---|---|
| SCREEN (Pan et al., 2024) | Structure, sequence embeddings, and evolutionary conservation | Graph Neural Network which employs contrastive learning on EC numbers |
| AEGAN (Shen et al., 2023) | Structure, Atchley Factors, PSSM conservation scores | Integrated graphical representations are generated using local protein networks and an amino acid spatial model with conservation scores and Atchley Factors included. A Graph Neural Network is trained to predict catalytic residues |
| PREvaIL (Song et al., 2018) | Structure and amino acid sequence | Generates a range of sequence- and structure-based features representing conservation, structural descriptions, and network properties of the enzyme, and then performs feature selection and trains Random Forest models for prediction |
| CRpred (Zhang et al., 2008) | Sequence | Generates and predicts sequence-based features, including PSSM conservation scores and trains an SVM for prediction |