Introduction

The diverse roles proteins play in nature necessitate their involvement in varied bio-functions, such as catalysis, binding, scaffolding, and transport (Frauenfelder and McMahon, 1998). Advances in science and technology have positioned proteins, particularly those with desired catalytic and binding properties, to pivotal positions in medicine, biology, and scientific research. While wild-type proteins exhibit optimal bio-functionality in their native environments, industrial applications often demand adaptation to conditions such as high temperature, high pressure, strong acidity, and strong alkalinity. The constrained efficacy of proteins in meeting the stringent requirements of industrial functioning environments hinders their widespread applications (Woodley, 2013; Ismail et al., 2021; Zheng et al., 2022; Jiang et al., 2023). Consequently, there arises a need to engineer these natural proteins to enhance their functionalities, aligning them with the demands of both industrial and scientific applications.

Analyzing the relationship between protein sequence and function yields valuable insights for engineering proteins with new or enhanced functions in synthetic biology. The intrinsic goal of protein engineering is to unveil highly functional sequences by navigating the intricate, high-dimensional surface of the fitness landscape (Yang et al., 2019), which delineates the relationship between amino acid (AA) sequences and the desired functional state. Predicting the effects of AA substitutions, insertions, and deletions (Riesselman et al., 2018; Gray et al., 2018; Notin et al., 2022) in proteins, i.e., mutation, thus allows researchers to dissect how changes in the AA sequence can impact the protein’s catalytic efficiency, stability, and binding affinity (Shin et al., 2021). However, the extensive, non-contiguous search space of AA combinations poses challenges for conventional methods, such as rational design (Aprile et al., 2020) or directed evolution (Wang et al., 2021), in efficiently and effectively identifying protein variants with desired properties. In response, deep learning has emerged as a promising solution that proposes favorable protein variants with high fitness.

Deep learning approaches have been instrumental in advancing scientific insights into proteins, predominantly categorized into sequence-based and structure-based methods. Autoregressive protein language models (Notin et al., 2022; Madani et al., 2023) interpret AA sequences as raw text to generate with self-attention mechanisms (Vaswani et al., 2017). Alternatively, masked language modeling objectives develop attention patterns that correspond to the residue-residue contact map of the protein (Meier et al., 2021; Rao et al., 2021a; Rives et al., 2021; Vig et al., 2021; Brandes et al., 2022; Lin et al., 2023). Other methods start from a multiple sequence alignment, summarizing the evolutionary patterns in target proteins (Riesselman et al., 2018; Frazer et al., 2021; Rao et al., 2021b). These methods result in a strong capacity for discovering the hidden protein space. However, the many-to-one relationship of both sequence-to-structure and structure-to-function projections requires excessive training input or substantial learning resources, which raises concerns regarding the efficiency of these pathways when navigating the complexities of the vast sequence space associated with a target function. Moreover, the overlooked local environment of proteins hinders the model’s ability to capture structure-sensitive properties that impact protein’s thermostability, interaction with substrates, and catalytic process (Zhao et al., 2022; Koehler Leman et al., 2023). Alternatively, structure-based modeling methods serve as an effective enhancement to complement sequence-oriented inferences for proteins with their local environment (Torng and Altman, 2019; Jones et al., 2021; Zhou et al., 2023b; Yi et al., 2023). Given that core mutations often induce functional defects through subtle disruptions to structure or dynamics (Roscoe et al., 2013), incorporating protein geometry into the learning process can offer valuable insights into stabilizing protein functioning. Recent efforts have been made to encode geometric information of proteins for topology-sensitive tasks such as molecule binding (Jin et al., 2021; Myung et al., 2022; Kong et al., 2023) and protein properties prediction (Zhang et al., 2022). Nevertheless, structure encoders fall short in capturing non-local connections for AAs beyond their contact region and overlook correlations that do not conform to the ‘structure-function determination’ heuristic.

There is a pressing need to develop a novel framework that overcomes limitations inherent in individual implementations of sequence or structure-based investigations. To this end, we introduce ProtSSN to assimilate the semantics and topology of Proteins from their Sequence and Structure with deep neural Networks. ProtSSN incorporates the intermediate state of protein structures and facilitates the discovery of an efficient and effective trajectory for mapping protein sequences to functionalities. The developed model extends the generalization and robustness of self-supervised protein language models while maintaining low computational costs, thereby facilitating self-supervised training and task-specific customization. A funnel-shaped learning pipeline, as depicted in Fig. 1, is designed due to the limited availability of protein crystal structures compared to observed protein sequences. Initially, the linguistic embedding establishes the semantic and grammatical rules in AA chains by inspecting millions of protein sequences. The topological embedding then enhances the interactions among locally connected AAs. Since a geometry can be placed or observed by different angles and positions in space, we represent proteins’ topology by graphs and enhance the model’s robustness and efficiency with a rotation and translation equivariant graph representation learning scheme.

Figure 1.

The pre-trained ProtSSN demonstrates its feasibility across a broad range of mutation effect prediction benchmarks covering catalysis, interaction, and thermostability. Mutation effects prediction serves as a common in silico assessment for evaluating the capability of an established deep learning framework in identifying favorable protein variants. With the continuous evolution of deep mutation scanning techniques and other high-throughput technologies, benchmark datasets have been curated to document fitness changes in mutants across diverse proteins with varying degrees of combinatorial and random modifications (Moal and Fernández-Recio, 2012; Nikam et al.,2021; Veleckỳ et al., 2022). ProteinGym v1 (Notin et al., 2023) comprehends fitness scoring of mutants to catalytic activity and binding affinity across over 200 well-studied proteins of varying deep mutational scanning (DMS) assays and taxa. Each protein includes tens of thousands of mutants documented from high-throughput screening techniques. While ProteinGym v1 initiates the most comprehensive benchmark dataset for mutants toward different properties enhancement, these fitness scores are normalized and simplified into several classes without further specifications. For instance, a good portion of protein assays is scored by their stability, which, in practice, is further categorized into thermostability, photostability, pH-stability, etc. Moreover, these stability properties should be discussed under certain environmental conditions, such as pH and temperature. Unfortunately, these detailed attributes are excluded in ProteinGym v1. Since a trade-off relationship emerges between activity and thermostability (Zheng et al., 2022), protein engineering in practice extends beyond enhancing catalytic activities to maintaining or increasing thermostability for prolonged functional lifetimes and applications under elevated temperatures and chemical conditions (Liu et al., 2019). Consequently, there is a need to enrich mutation effect prediction benchmarks that evaluate the model’s efficacy in capturing variant effects concerning thermostability under distinct experimental conditions. In this study, we address the mutation effect prediction tasks on thermostability with DTm and DDG, two new single-site mutation benchmarks that measure thermostability using ΔTm and ΔΔG values, respectively. Both benchmarks group experimental assays based on the protein-condition combinations. These two datasets supplement the publicly available DMS assay benchmarks and facilitate ready-to-use assessment for future deep learning methods toward protein thermostability enhancement.

This work fulfills practical necessities in the development of deep learning methods for protein engineering from three distinct perspectives. Firstly, the developed ProtSSN employs self-supervised learning during training to obviate the necessity for additional supervision in down-stream tasks. This zero-shot scenario is desirable due to the scarcity of experimental results, as well as the ‘cold-start’ situation common in many wet lab experiments. Secondly, the trained model furnishes robust and meaningful approximations to the joint distribution of the entire amino acid chain, thereby augmenting the epistatic effect (Sarkisyan et al., 2016; Khersonsky et al., 2018) in deep mutations. This augmentation stems from the model’s consideration of the nonlinear combinatorial effects of amino acid sites. Additionally, we complement the existing dataset of DMS assay with additional benchmark assays related to protein-environment interactions, specifically focusing on thermostability–an integral objective in numerous protein engineering projects.

Materials and Methods

This section starts by introducing the three mutational scanning benchmark datasets used in this study, including an open benchmark dataset ProteinGym v1 and two new datasets on protein ther-mostability, i.e., DTm and DDG. Section establishes the proposed ProtSSN for protein sequences and structures encoding. The overall pipeline of model training and inference is presented in Section . The experimental setup for model evaluation is provided in Section .

Dataset

To evaluate model performance in predicting mutation effects, we compare ProtSSN with SOTA baselines on ProteinGym-v1, the largest open benchmark for DMS assays. We also introduce two novel benchmarks, DTm and DDG, for assessing protein thermostability under consistent control environments. Details of the new datasets are provided in Table 1. Considering the significant influence of experimental conditions on protein temperature, we explicitly note the pH environment in each assay. Refer to the following paragraphs for more details.

Table 1.

ProteinGym

The assessment of ProtSSN for deep mutation effects prediction is conducted on ProteinGym v1 (Notin et al., 2023). It is the most extensive protein substitution benchmark comprising 217 assays and more than 2.5M mutants. These DMS assays cover a broad spectrum of functional properties (e.g., ligand binding, aggregation, viral replication, and drug resistance) and span various protein families (e.g., kinases, ion channel proteins, transcription factors, and tumor suppressors) across different taxa (e.g., humans, other eukaryotes, prokaryotes, and viruses).

DTm

The assays in the first novel thermostability benchmark are sourced from single-site mutations in ProThermDB (Nikam et al., 2021). Each assay is named in the format ‘UniProt_ID-pH’. For instance, ‘O00095-8.0’ signifies mutations conducted and evaluated under pH 8.0 for protein O00095. The attributes include ‘mutant’, ‘score’, and ‘UniProt_ID’, with at least 10 mutants to ensure a meaningful evaluation. To concentrate on single-site mutations compared to wild-type proteins, we exclude records with continuous mutations. To avoid dealing with partial or misaligned protein structures resulting from incomplete wet experimental procedures, we employ UniProt ID to pair protein sequences with folded structures predicted by AlphaFold2 (Jumper et al., 2021) from https://alphafold.ebi.ac.uk. In total, DTm consists of 66 such protein-environment pairs and 2, 520 mutations.

DDG

The second novel thermostability benchmark is sourced from DDMut (Zhou et al., 2023c), where mutants are paired with their PDB documents. We thus employ crystal structures for proteins in DDG and conduct additional preprocessing steps for environmental consistency. In particular, we group mutants of the same protein under identical experimental conditions and examine the alignment between mutation sites and sequences extracted from crystal structures. We discard mutant records that cannot be aligned with the sequence loaded from the associated PDB document. After removing data points with fewer than 10 mutations, DDG comprises 36 data points and 829 mutations. As before, assays are named in the format ‘PDB_ID-pH-temperature’ to indicate the chemical condition (pH) and temperature (in Celsius) of the experiment on the protein. An example assay content is provided in Fig. 2.

Figure 2.

Proposed Method

Labeled data are usually scarce in biomolecule research, which demands designing a general self-supervised model for predicting variant effects on unknown proteins and protein functions. We propose ProtSSN with a pre-training scheme that recovers AA types from a given protein backbone structure with noisy local environment observations.

Multi-level Protein Representation

Protein Primary Structure (Noisy)

Denote a given observation with AA type . We assume this arbitrary observation is undergoing random mutation toward a stable state. The training objective is to learn a revised state v that is less likely to be eliminated by natural selection due to unfavorable properties such as instability or inability to fold. The perturbed observation is defined by a multinomial distribution.

where an AA in a protein chain has a chance of p to mutate to one of 20 AAs (including itself) following a replacement distribution Θ(·) and (1 − p) of remaining unchanged. We consider p as a tunable parameter and define Θ(·) by the frequency of AA types observed in wild-type proteins.

Protein Tertiary Structure

The geometry of a protein is described by 𝒢= (𝒱,ε?, W _V, W _E, X_V ), which is a residue graph based on the k-nearest neighbor (kNN). Each node v_i ∈ 𝒱 represents an AA in the protein connected to up to k closest nodes within 30Å (Zhou et al., 2023a). Node attributes W _V are hidden semantic embeddings of AA types extracted by the semantic encoder (see Section ), and edge attributes W _E ∈ ℝ⁹³ feature relationships of connected nodes based on inter-atomic distances, local N-C positions, and sequential position encoding. Additionally, X_V records 3D coordinates of AAs in the Euclidean space, which plays a crucial role in the subsequent geometric embedding stage to preserve roto-translation equivariance.

Semantic Encoding of Global AA Contacts

Although it is generally believed by the community that a protein’s sequence determines its biological function via the folded structure, following strictly to this singular pathway risks overlooking other unobserved yet potentially influential inter-atomic communications impacting protein fitness. Our proposed ProtSSN thus includes pairwise relationships for residues through an analysis of proteins’ primary structure from and embeds them into hidden representations W _V for residues. At each update, information and representations of the noisy AA sequence are encoded from the noisy input via an Evolutionary Scale Modeling method, ESM-2 (Lin et al., 2023), which employs a BERT-style masked language modeling objective. This objective predicts the identity of randomly selected AAs in a protein sequence by observing their context within the remainder of the sequence. Alternative semantic encoding strategies are discussed in Section .

Geometrical Encoding of Local AA Contacts

Proteins are structured in 3D space, which requires the geometric encoder to possess roto-translation equivariance to node positions as well as permutation invariant to node attributes. This design is vital to avoid the implementation of costly data augmentation strategies. We practice Equivariant Graph Neural Networks (EGNN) (Satorras et al., 2021) to acquire the hidden node representation and node coordinates at the l + 1th layer

In these equations, represents the input edge attribute on 𝒱_ij, which is not updated by the network. The propagation rules ϕ_e, ϕ_x and ϕ_v are defined by differentiable functions, e.g.multi-layer perceptrons (MLPs). The final hidden representation on nodes embeds the microenvironment and local topology of AAs, and it will be carried on by readout layers for label predictions.

Model Pipeline

ProtSSN is designed for protein engineering with a self-supervised learning scheme. The model is capable of conducting zero-shot variant effect prediction on an unknown protein, and it can generate the joint distribution for the type of all AAs along the protein sequence of interest. This process accounts for the epistatic effect and concurrently returns all AA sites in a sequence. Below, we detail the workflow for training the zero-shot model and scoring the effect of a specific mutant.

Training

The fundamental model architecture cascades a frozen sequence encoding module and a trainable tertiary structure encoder. The protein is represented by its sequence and structure characteristics, where the former is treated as plain text and the latter is formalized as a kNN graph with model-encoded node attributes, handcrafted edge attributes, and 3D positions of the corresponding AAs. The node attributes are encoded by a protein language model, which analyses pairwise hidden relationships among AAs from the input protein sequence and produces a vector representation for each AA, i.e.,

with be the perturbed initial AA-type encoding. The language model LM_frozen(·), ESM-2 (Lin et al., 2023) for instance, has been pre-trained on a massive protein sequence database (e.g., UniRef50 (Suzek et al., 2015)) to understand the semantic and grammatical rules of wild-type proteins with high-dimensional AA-level long-short-range interactions. The independent perturbation on AA types follows 1 operates in every epoch, which requires constant updates on the sequence embedding. To further enhance the interaction of locally connected AAs in the node representation, a stack of L EGNN layers is implemented on the input graph 𝒢 to yield

During the pre-training phase for protein sequence recovery, the output layer ϕ(·) provides the probability of AAs in one of 20 types, i.e.,

for n AAs of the input protein. The model’s learnable parameters are refined by minimizing the cross-entropy of the recovered AAs with respect to the ground-truth AAs in wild-type proteins.

Inferencing

For a given mutant, its fitness score is derived from the joint distribution of the altered AA types on associated nodes that provides a preliminary assessment based on natural observations. We consider the AA type in the wild-type protein as a reference state and compare it with the predicted probability of AAs at the mutated site. Formally, for a mutant with mutated sites 𝒯 (|𝒯 |≥ 1), we define its fitness score using the corresponding log-odds-ratio, i.e.,

where y_t and v_t denote the mutated and wild-type type of the tth AA.

Experimental Protocol

We validate the efficacy of ProtSSN on zero-shot mutation effect prediction tasks. The performance is compared with other SOTA models of varying numbers of trainable parameters. The implementations are programmed with PyTorch-Geometric (ver 2.2.0) and PyTorch (ver 1.12.1) and executed on an NVIDIA^® Tesla A100 GPU with 6, 912 CUDA cores and 80GB HBM2 installed on an HPC cluster.

Training Setup

ProtSSN is pre-trained on a non-redundant subset of CATH v4.3.0 (Orengo et al., 1997) domains, which contains 30, 948 experimental protein structures with less than 40% sequence identity. We remove proteins that exceed 2, 000 AAs in length for efficiency. For each kNN protein graph, node features are extracted and updated by a frozen ESM2-t33 prefix model (Lin et al., 2023). Protein topology is inferred by a 6-layer EGNN (Satorras et al., 2021) with the hidden dimension tuned from {512, 768, 1280}. Adam (Kingma and Ba, 2015) is used for backpropagation with the learning rate set to 0.0001. To avoid training instability or CUDA out-of-memory, we limit the maximum input to 8, 192 AA tokens per batch, constituting approximately 32 graphs.

Baseline Methods

We undertake an extensive comparison with baseline methods of self-supervised sequence or structure-based models on the fitness of mutation effects prediction (Table 2). Sequence models employ position embedding strategies such as autoregression (RITA (Hesslow et al., 2022), Tranception (Notin et al., 2022), and ProGen2 (Nijkamp et al., 2023)), masked language modeling (ESM-1b (Rives et al., 2021), ESM-1v (Meier et al., 2021), and ESM2 (Lin et al., 2023)), a combination of the both (ProtTrans (Elnaggar et al., 2021)), and convolutional modeling (CARP (Yang et al., 2022)). Additionally, we also compare with ESM-if1 (Hsu et al., 2022) which incorporates masked language modeling objectives with GVP (Jing et al., 2020), as well as MIF-ST (Yang et al., 2023) with both sequence and structure encoding.

Table 2.

Benchmark Datasets

We conduct a comprehensive comparison of deep learning predictions on hundreds of DMS assays concerning different protein types and biochemical assays. We prioritize experimentally measured properties that possess a monotonic relationship with protein fitness, such as catalytic activity, binding affinity, expression, and thermostability. In particular, ProteinGym v1 (Notin et al., 2023) groups 5 sets of protein properties from 217 assays. Two new datasets named DTm and DDG examine 102 environment-specific assays with thermostability scores. In total, 319 assays with around 2.5 million mutants are scored. To better understand the novelty of the proposed ProtSSN, we designed additional investigations on ProteinGym v0 (Notin et al., 2022), an early-release version of ProteinGym with 1.5M missense variants across 86 assays (excluding one assay that failed to be folded by both AlphaFold2 and ESMFold, i.e., A0A140D2T1_ZIKV_Sourisseau_growth_2019 ).

Evaluation Metric

We evaluate the performance of pre-trained models on a diverse set of proteins and protein functions using Spearman’s ρ correlation that measures the strength and direction of the monotonic relationship between two ranked sequences, i.e., experimentally evaluated mutants and model-inferred mutants (Meier et al., 2021; Notin et al., 2022; Frazer et al., 2021; Zhou et al., 2023a). This non-parametric rank measure is robust to outliers and asymmetry in mutation scores, and it does not assume any specific distribution of mutation fitness scores. The scale of ρ ranges from −1 to 1, indicating the negative or positive correlation of the predicted sequence to the ground truth. The prediction is preferred if its ρ to the experimentally examined ground truth is close to 1.

Results

Variant Effect Prediction

We commence our investigation by assessing the predictive performance of ProtSSN on four bench-mark datasets using state-of-the-art (SOTA) protein learning models. We deploy different versions of ProtSSN that learns kNN graphs with k ∈ {10, 20, 30} and h ∈ {512, 768, 1, 280} hidden neurons in each of the 6 EGNN layers. For all baselines, their performance on DTm and DDG are reproduced with the official implementation listed in Table 2, and the scores on ProteinGym v1 are retrieved from https://proteingym.org/benchmarks by Notin et al. (2023). As visualized in Fig. 3 and reported in Table 3, ProtSSN demonstrated exceptional predictive performance with a significantly smaller number of trainable parameters when predicting the function and thermostability of mutants.

Table 3.

Figure 3.

Compared to protein language models (Fig. 3, colored circles), ProtSSN benefits from abundant structure information to more accurately capture the overall protein characteristics, resulting in enhanced thermostability and thus achieving higher correlation scores on DTm and DDG. This is attributed to the compactness of the overall architecture as well as the presence of stable local structures such as alpha helices and beta sheets, both of which are crucial to protein thermostability by providing a framework that resists unfolding at elevated temperatures (Robinson-Rechavi et al., 2006). Consequently, the other two structures-involved models, i.e., ESM-if1 and MIF-ST, also exhibit higher performance on the two thermostability benchmarks.

While sequence-based analysis cannot entirely replace the role of structure-based analysis, compared to a fully structure-based deep learning method, a protein language model is more likely to capture sufficient information from sequences by increasing the model scale, i.e., the number of trainable parameters. This observation aligns with the mechanism of protein functionality, where local structures (such as binding pockets and catalytic pockets) and the overall structure (such as conformational changes and allosteric effect) are both crucial for binding and catalysis (Sheng et al., 2014). Unlike the thermostability benchmarks, the discrepancies between structure-involved models and protein language models are mitigated in ProteinGym v1 by increasing the model scale. Overall, structure-based and sequence-based methods are suitable for different types of assays, but a comprehensive framework, such as ProtSSN, exhibits the best performance.

Ablation Study

The effectiveness of ProtSSN with different modular designs is examined by its performance on the early released version ProteinGym v0 (Notin et al., 2022) with 86 DMS assays. The complete comparison of baseline comparison on ProteinGym v0 is detailed in Table S1 in SI. For the ablation study, the results are summarized in Fig. 4(a)-(d) and detailed in Table S2, where each record is subject to a combination of key arguments under investigation. For instance, in the top orange box of Fig. 4(a), we report all ablation results that utilize 6 EGNN layers for graph convolution, regardless of the different scales of ESM-2 or the definitions of node attributes. For all modules investigated in this section, we separately discuss their influence on predicting mutation effects when modifying a single site or an arbitrary number of sites. The two cases are marked on the y-axis by single’ and ‘all’, respectively.

Figure 4.

Inclusion of Roto-Translation Equivariance

We assess the effect of incorporating rotation and translation equivariance for protein geometric and topological encoding. Three types of graph convolutions are compared, including GCN (Kipf and Welling, 2017), GAT (Veličković et al., 2018), and EGNN (Satorras et al., 2021). The first two are classic non-equivariant graph convolutional methods, and the last one, which we apply in the main algorithm, preserves roto-translation equivariance. We fix the number of EGNN layers to 6 and examine the performance of the other two methods with either 4 or 6 layers. We find that integrating equivariance when embedding protein geometry significantly.

Sequence Encoding

We next investigate the benefits of defining data-driven node attributes for protein representation learning. We compare the performance of models trained on two sets of graph inputs: the first set defines its AA node attributes through trained ESM2 (Lin et al., 2023), while the second set uses one-hot encoded AA types for each node. A clear advantage of using hidden representations by prefix models over hardcoded attributes is evident from the results presented in Fig. 4(b).

Depth of EGNN

Although graph neural networks can extract topological information from geometric inputs, it is vital to select an appropriate number of layers for the module to deliver the most expressive node representation without encountering the oversmoothing problem. We investigate a wide range of choices for the number of EGNN layers among {6, 12, 18}. As reported in Fig. 4(c), embedding graph topology with deeper networks does not lead to performance improvements. A moderate choice of 6 EGNN layers is sufficient for our learning task.

Scale of ESM

We also evaluate ProtSSN on different choices of language embedding dimensions to study the trade-off between computational cost and input richness. Various scales of prefix models, including {8, 150, 650, 3000} millions of parameters, have been applied to produce different sequential embeddings with {320, 640, 1280, 2560} dimensions, respectively. Fig. 4(d) reveals a clear preference for ESM2-t33, which employs 650 million parameters to achieve optimal model performance with the best stability. Notably, a higher dimension and richer semantic expression do not always yield better performance. A performance degradation is observed at ESM2 with 3 billion parameters.

AA Perturbation Strategies

ProtSSN utilizes noise sampling at each epoch on the node attributes to emulate random mutations in nature. This introduction of noise directly affects the node attribute input to the graphs. Alongside the ‘mutate-then-recode’ method we implemented in the main algorithm, we examined four additional strategies to perturb the input data during training. The construction of these ad-ditional strategies is detailed below, and their corresponding Spearman’s ρ on ProteinGym v0 is in Fig. 4(e).

Global Average

Suppose the encoded sequential representation of a node is predominantly determined by its AA type. In essence, the protein sequence encoder will return similar embeddings for nodes of the same AA type, albeit at different positions within a protein. With this premise, the first strategy replaces the node embedding for perturbed nodes with the average representations of the same AA types. For example, a random node v_i in a protein is Glycine. If it is altered in the current epoch, is designated as the average sequential embedding of all other Glycines.

Sliding Window

The presumption in the previous method neither aligns with the algorithmic design nor biological heuristics. Self-attention discerns the interaction of the central token with its neighbors across the entire document (protein). The representation of an AA is influenced by both its neighbors’ types and locations. Thus, averaging embeddings of AAs from varying positions is likely to forfeit positional information of the modified AA. For this purpose, we consider a sliding window of size 3 along the protein sequence.

Gaussian Noise

The third option regards node embeddings of AA as hidden vectors and imposes white noise on the vector values. We define the mean= 0 and variance= 0.5 on the noise, making the revised node representation .

Mask

Finally, we employ the masking technique prevalent in masked language modeling and substitute the perturbed AA token with a special <mask> token. The prefix language model will interpret it as a novel token and employ self-attention mechanisms to assign a new representation to it. We utilize the same hyperparameter settings as that of BERT (Devlin et al., 2018) and choose 15% of tokens per iteration as potentially influenced subjects. Specifically, 80% of these subjects are replaced with <mask>, 10% of them are replaced randomly (to other AA types), and the remaining 10% stay unchanged.

Folding Methods

The analysis of protein structure through topological embeddings poses challenges due to the inherent difficulty in obtaining accurate structural observations through experimental techniques such as NMR (Wüthrich, 1990) and Cryo-EM (Elmlund et al., 2017). As a substitution, we use in silico folding models such as AlphaFold2 (Jumper et al., 2021) and ESMFold (Lin et al., 2023). Although these folding methods may introduce additional errors due to the inherent uncertainty or inaccuracy in structural predictions, we investigate the robustness of our model in the face of potential inaccuracies introduced by these folding strategies.

Table 4 examines the impact of different folding strategies on the performance of DMS predictions for ESM-if1, MIF-ST, and ProtSSN. Although the performance based on ESMFold-derived protein structures generally lags behind structures folded by AlphaFold2, the minimal difference between the two strategies across all four benchmarks validates the robustness of our ProtSSN in response to variations in input protein structures. We deliberately divide the assays in ProteinGym v0 (Notin et al., 2022) into two groups of interaction (e.g., binding and stability) and catalysis (e.g., activity), as the former is believed more sensitive to the structure of the protein (Robertson and Murphy, 1997).

Table 4.

ProtSSN emerges as the most robust method, maintaining consistent performance across different folding strategies, which underscores the importance of our model’s resilience to variations in input protein structures. The reported performance difference for both ESM-if1 and MIF-ST is 3−5 times higher than that of ProtSSN. The inconsistency between the optimal results for DDG and the scores reported in Table 3 comes from the utilization of crystal structures in the main results. Another noteworthy observation pertains to the performance of ESM-if1 and MIF-ST on DDG. In this case, predictions based on ESMFold surpass those based on AlphaFold2, even outperforming predictions derived from crystal structures. However, the superior performance of these methods with ESMFold hinges on inaccurate protein structures, rendering them less reliable.

Performance Enhancement with MSA and Ensemble

We extend the evaluation on ProteinGym v0 to include a comparison of our zero-shot ProtSSN with few-shot learning methods that leverage MSA information of proteins (Meier et al., 2021). The details are reported in Table 5, where the performance of baseline methods are retrieved from https://github.com/OATML-Markslab/ProteinGym. For ProtSSN, we employ an ensemble of nine models with varying sizes of kNN graphs (k ∈ {10, 20, 30}) and EGNN hidden neurons ({512, 768, 1280}), as discussed previously. The MSA depth reflects the number of MSA sequences that can be found for the protein, which influences the quality of MSA-based methods.

Table 5.

Table 6.

ProtSSN demonstrates superior performance among non-MSA zero-shot methods in predicting both single-site and multi-site mutants, underscoring its pivotal role in guiding protein engineering towards deep mutation. Furthermore, ProtSSN outperforms MSA-based methods across taxa, except for viruses, where precise identification of conserved sites is crucial for positive mutations, especially in spike proteins (Katoh and Standley, 2021; Mercurio et al., 2021). In essence, ProtSSN provides an efficient and effective solution for variant effects prediction when compared to both non-MSA and MSA-based methods. Note that although MSA-based methods generally consume fewer trainable parameters than non-MSA methods, they incur significantly higher costs in terms of search time and memory usage. Moreover, the inference performance of MSA-based methods relies heavily on the quality of the input MSA, where the additionally involved variance makes impacts the stability of the model performance.

Discussion and Conclusion

The development of modern computational methodologies for protein engineering is a crucial facet of in silico protein design. Effectively assessing the fitness of protein mutants not only supports cost-efficient experimental validations but also guides the modification of proteins toward enhanced or novel functions. Recent deep learning solutions employ a common strategy that involves establishing a hidden protein representation and masking potential mutation sites to recommend plausible amino acids. Previous research has primarily focused on extracting protein representations from either their sequential or structural modalities, with many treating the prediction of mutational effects merely as a secondary task following inverse folding or de novo protein design. These approaches often overlook the importance of comprehensive investigation on both global and local perspectives of amino acid interaction that are critical for featuring protein functions. Furthermore, these methods hardly tailored model designs for suggesting mutations, despite the significance of this type of tasks. In this work, we introduce ProtSSN, a denoising framework that effectively cascades protein sequence and structure embedding for predicting mutational effects. Both the protein language model and equivariant graph neural network are employed to encode the global interaction of amino acids with geometry-aware local enhancements. We have extensively validated the efficacy of ProtSSN across various protein function assays and taxa, including two thermostability datasets prepared by ourselves. Our approach consistently demonstrates substantial promise for protein engineering applications, such as facilitating the design of mutation sequences with enhanced interaction and enzymatic activity.