Method for predicting protein mutation stability change based on integrated convolutional neural network model and regression hierarchical training

Abstract

The protein mutation stability change prediction method based on integrated convolutional neural network model and regression hierarchical training belongs to the technical field of protein stability change prediction. First, a training set and a test set are constructed, and four data sets are collected and arranged. Then the following four steps are carried out: first, the data is enhanced based on antisymmetry, and the sample data of stable and unstable mutations in the training set is balanced; second, protein features are extracted; third, a regression hierarchical sampling strategy is used to train the model; fourth, an integrated model based on multiple CNN sub-models is constructed for protein mutation stability change prediction. The present application first introduces spatial neighborhood evolutionary information; for the first time, a regression hierarchical sampling strategy is proposed and successfully applied in the training process of the model, effectively improving the prediction ability of the model for extreme ΔΔG; the trained CNN sub-models are combined, and the information of all samples in the training set is fully utilized, which is conducive to enhancing the generalization ability of the model.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of protein stability change prediction, and is a protein mutation stability change prediction method based on protein sequence and structure information. BACKGROUND

[0002] Many human genetic diseases are directly related to missense mutations in proteins [1] . Missense mutations usually lead to an increase (stabilizing mutations) or decrease (destabilizing mutations) in protein stability, thereby affecting its function and triggering diseases [2] . Currently, accurately predicting the stability change caused by protein mutations is an important and challenging task. Generally, the stability change of protein caused by mutation is quantified by measuring the change in folding free energy (ΔΔG) before and after mutation [1] . Currently, experimental methods such as differential scanning calorimetry and thermal fluorescence spectroscopy can detect protein stability changes [3] , but these methods are time-consuming and tedious. Therefore, there is an urgent need to develop an efficient computational method to predict the stability change of protein caused by mutation.

[0003] Several theoretical methods have been developed to predict protein stability changes (ΔΔG), which can be mainly divided into three categories: energy-based, knowledge-based and machine learning-based methods [4] . Free energy perturbation (FEP) [5] , an early developed method based on molecular dynamics simulation, is generally considered to have high accuracy. FoldX model is a traditional energy-based method that uses empirical force fields and energy terms to calculate the effect of mutation on protein stability [6] . Although these two methods have relatively good performance, they are computationally complex and time-consuming. For knowledge-based methods, there are two well-known tools, SDM [7] and DdGUN [8] . SDM first used conformationally constrained environment-specific substitution tables (ESSTs) to calculate the change in protein stability before and after mutation, while DdGUN only used evolutionary information to calculate the difference in stability. Both tools use statistical methods, which have limited generalization ability. In recent years, many machine learning-based prediction methods have been developed, including MAESTRO [9] , mCSM

[10] , DUET

[11] , I-Mutant

[12] , iStable

[13] and ThermoNet

[14] and so on. These methods use traditional tree models or deep neural networks and consider various features, including sequence, structure, and physicochemical features, etc., to improve the generalization ability of the model to some extent.

[0004] Most of the methods for predicting the stability change of protein mutations have two common problems: one is poor prediction ability for mutations with stability mutations; the second is poor prediction ability for mutations with extreme ΔΔG values (highly unstable or highly stable mutations). For the first problem, the main reason is that the training set used by the model is unbalanced (more than 75% of the mutations in the existing data set are unstable

[15] ), which will cause the prediction of the model to be biased. At present, most prediction methods ignore this problem. The stability change (ΔΔG) of protein caused by mutation has an important property called anti-symmetry, that is, the Gibbs free energy change of amino acid A to amino acid B (direct mutation) at a certain position of the protein is A→B = ΔG A - ΔG B , and the Gibbs free energy change of amino acid B to amino acid A (reverse mutation) is B→A = ΔG B - ΔG A , so we can get A→B = - ΔΔG B→A

[16] . We consider using this property to construct a balanced training set, so as to improve the prediction bias caused by unbalanced data. For the second problem, the main reason is that generally there are few samples of mutations with extreme ΔΔG values, so the information that the model can learn is limited. How to make full use of the information contained in these few samples in the training process is very important for the model to accurately predict extreme values. To our knowledge, existing methods ignore this point, so there is no measure to solve this problem. It is worth noting that generally mutations with extreme ΔΔG values are more likely to have a greater impact on the structure of the protein and thus cause diseases. In order to solve this problem, we design a regression stratified sampling training method, which will help to improve the prediction ability of the model for extreme ΔΔG values.

[0005] In addition, feature extraction is also crucial for the model. We considered physicochemical properties, Position Specific Scoring Matrix (PSSM), secondary structure features, as well as the depth index and protrusion index. Moreover, we recently proposed a new evolutionary information encoding method, Spacial Neighbor-Based Position-Specific Scoring Matrix (SNB-PSSM), which takes into account the characteristics of spatially co-evolving residues (averaging over spatially neighboring residues) when encoding the evolutionary information of residues, and achieved good results in protein-RNA binding site prediction.

[17] In addition to the above traditional features, we believe that further consideration of topological information can enhance the performance of the model. Persistent homology, as a branch of algebraic topology, provides an effective method to describe the topological structure of proteins

[18] , and has been successfully applied to many aspects, such as protein-ligand binding affinity prediction, protein-protein interaction prediction, and drug virtual screening [19-21] . Therefore, we believe that topological features based on persistent homology will also help improve the model's ability to predict protein stability changes.

[0006] In recent years, deep learning has been widely applied in the field of protein structure and function prediction, such as protein structure prediction, functional annotation, interaction prediction, and drug design. Among them, Convolutional Neural Networks (CNN) is more efficient and accurate in processing large-scale protein data, and has the advantages of translation invariance and parameter sharing, which performs well in capturing local features in data. Most deep learning will find the optimal model in the training process as the final prediction model, but the samples in the validation set do not participate in the training of the model during training, resulting in the loss of information contained in this part of the sample. Ensemble Method

[22] can solve this problem well. Therefore, we consider adopting Ensemble Method to construct the CNN model framework.

[0007] In summary, we consider the antisymmetry to solve the problem of data imbalance; by extracting protein sequence and structure features, and newly introducing spatial neighbor evolutionary information and topological features based on persistent homology, we fully consider the features related to protein stability; we propose a stratification sampling strategy called Regression Stratification cross-validation (RScv) to solve the problem of poor prediction effect for mutations with extreme ΔΔG values; we propose an ensemble method to combine the trained CNN sub-models to fully utilize the information of all samples in the training set. Finally, we construct an unbiased model for predicting the stability change of protein mutations based on the ensemble convolutional neural network model and regression stratification training. This method can provide important information for protein engineering and related drug design.

[0008] REFERENCES

[0009] [1] STEFL S, NISHI H, PETUKH M, et al. Molecular Mechanisms of Disease-Causing Missense Mutations [J]. Journal of Molecular Biology, 2013, 425(21): 3919-3936.

[0010] [2] BANERJEE A, MITRA P. Estimating the Effect of Single-Point Mutations on Protein Thermodynamic Stability and Analyzing the Mutation Landscape of the p53 Protein [J]. J Chem Inf Model, 2020, 60(6): 3315-3323.

[0011] [3] SANAVIA T, BIROLO G, MONTANUCCI L, et al. Limitations and challenges in protein stability prediction upon genome variations: towards future applications in precision medicine [J]. Comput Struct Biotechnol J, 2020, 18: 1968-1979.

[0012] [4] PAN Q, NGUYEN T B, ASCHER D B, et al. Systematic evaluation of computational tools to predict the effects of mutations on protein stability in the absence of experimental structures [J]. Brief Bioinform, 2022, 23(2): bbae025.

[0013] [5] ZWANZIG R. High-temperature equation of state by a perturbation method [J]. The Journal of Chemical Physics, 1954, 8(22): 1420-1426.

[0014] [6] GUEROIS R, NIELSEN J E, SERRANO L. Predicting changes in the stability of proteins and protein complexes: a study of more than 1000 mutations [J]. J Mol Biol, 2002, 320(2): 369-387.

[0015] [7] PANDURANGAN A P, OCHOA-MONTANO B, ASCHER D B, et al. SDM: a server for predicting effects of mutations on protein stability [J]. Nucleic Acids Res, 2017, 45(W1): W229-W235.

[0016] [8] MONTANUCCI L, CAPRIOTTI E, FRANK Y, et al. DDGun: an untrained method for the prediction of protein stability changes upon single and multiple point variations [J]. BMC Bioinformatics, 2019, 20(Suppl 14): 335.

[0017] [9] LAIMER J, HOFER H, FRITZ M, et al. MAESTRO--multi agent stability prediction upon point mutations [J]. BMC Bioinformatics, 2015, 16: 116.

[0018]

[10] PIRESD E, ASCHER D B, BLUNDELL T L. mCSM: predicting the effects of mutations in proteins using graph-based signatures [J]. Bioinformatics, 2014, 30(3): 335-342.

[0019]

[11] PIRESD E, ASCHER D B, BLUNDELL T L. DUET: a server for predicting effects of mutations on protein stability using an integrated computational approach [J]. Nucleic Acids Res, 2014, 42 (Web Server issue): W314-W319.

[0020]

[12] CAPRIOTTI E, FARISELLI P, CASADIO R. I-Mutant2.0: predicting stability changes upon mutation from the protein sequence or structure [J]. Nucleic Acids Res, 2005, 33 (Web Server issue): W306-W310.

[0021]

[13] CHEN C W, LIN J, CHU Y W. iStable: off-the-shelf predictor integration for predicting protein stability changes [J]. BMC Bioinformatics, 2013, 14 Suppl 2: S5.

[0022]

[14] LI B, YANG Y T, CAPRA J A, et al. Predicting changes in protein thermodynamic stability upon point mutation with deep 3D convolutional neural networks [J]. PLoS Comput Biol, 2020, 16(11): e1008291.

[0023]

[15] TOKURIKI N, TAWFIK D S. Stability effects of mutations and protein evolvability [J]. Curr Opin Struct Biol, 2009, 19(5): 596-604.

[0024]

[16] THILTGEN G, GOLDSTEIN R A. Assessing predictors of changes in protein stability upon mutation using self-consistency [J]. PLoS One, 2012, 7(10): e46084.

[0025]

[17] LIU Y, GONG W, YANG Z, et al. SNB-PSSM: A spatial neighbor-based PSSM used for protein-RNA binding site prediction [J]. J Mol Recognit, 2021, 34(6): e2887.

[0026]

[18] ZOMORODIAN A, CARLSSON G. Computing persistent homology [J]. DISCRETE & COMPUTATIONAL GEOMETRY, 2005, 33(2): 249-274.

[0027]

[19] CANG Z, WEI G W. Integration of element specific persistent homology and machine learning for protein-ligand binding affinity prediction [J]. Int J Numer Method Biomed Eng, 2018, 34(2): e2914.

[0028]

[20] WANG M, CANG Z, WEI G W. A topology-based network tree for the prediction of protein-protein binding affinity changes following mutation [J]. Nat Mach Intell, 2020, 2(2): 116-123.

[0029]

[21] XIA K, WEI G W. Persistent homology analysis of protein structure, flexibility, and folding [J]. Int J Numer Method Biomed Eng, 2014, 30(8): 814-844.

[0030]

[22] DIETTERICH T G. Ensemble Methods in Machine Learning [C]. In: Multiple Classifier Systems, Berlin, Heidelberg, 2000: 1-15. SUMMARY

[0031] The purpose of the present application is to predict the stability change (ΔΔG) of a certain site amino acid mutation in a known protein, which helps to study the structure-function relationship of the protein and provides important information for protein engineering and drug design. The features of the method are: ① Based on the anti-symmetry, the data is enhanced and a balanced data set (containing the same number of stable and unstable mutations) is finally generated, which helps to solve the prediction bias problem of unstable mutation samples and can improve the ability to predict stable mutations. ② In terms of feature extraction, in addition to traditional protein sequence and structure features, a new evolutionary information encoding method, i.e., position-specific scoring matrix based on spatial neighbors, is first introduced, which is developed by our group (J Mol Recognit.2021 Jun 1;6(34):e2887). In addition, topological features based on persistent homology are also used to predict protein stability changes, which can play an important role in predicting protein stability. ③ Regression stratification sampling strategy is first proposed and applied to the prediction of protein mutation stability changes, which effectively improves the prediction ability of the model for extreme ΔΔG. ④ Unlike previous methods that use deep learning to predict protein mutation stability changes, we propose a new effective ensemble model containing multiple CNN sub-models based on ensemble method, which is helpful to enhance the generalization ability of the model.

[0032] The present application is based on an ensemble convolutional neural network model and a regression stratification cross-validation strategy, which is a protein mutation stability change prediction method (an ensemble convolutional neural network model with Regression Stratification cross-validation for Protein Stability changes upon Mutations Prediction, PMSPcnn) considering anti-symmetry and using spatial neighbor evolutionary information and persistent homology to extract topological features, including four steps: ① Based on anti-symmetry, the data is enhanced and the sample data of stable and unstable mutations in the training set is balanced, ② Protein features are extracted, ③ Regression stratification sampling strategy is used to train the model, and ④ an ensemble model based on multiple CNN sub-models is constructed for protein mutation stability change prediction (as shown in FIG. 2). Figure 1 )。

[0033] Before building the model, first construct the training set and test set:

[0034] Experimental data about protein mutation stability change (ΔΔG) from recently published literatures were collected: Q4321 dataset (Bioinformatics. 2016 Oct 1; 19(32): 2936-2946); S sym dataset (Bioinformatics. 2018 Nov 1; 21(34): 3659-3665); P53 dataset (Oncogene. 2000 Mar 2; 10(19): 1245-1256); Myoglobin dataset (J Exp Biol. 2018 Nov 1; 20(207): 3441-3446). Q4321 was used as training set, and the other three datasets S sym , P53 and Myoglobin were used as test sets. The samples in the training set that were also in the test set were removed, and then the proteins in the training set with E-value (an index to evaluate the similarity of two sequences, the smaller the E-value, the higher the similarity of the two sequences) less than 0.001 to the proteins in the test set were removed using BLAST method (https: / / ftp.ncbi.nlm.nih.gov / blast / executables / blast+ / LATEST / ). After processing, the training set (named as Q3211) finally contained 147 proteins, a total of 3211 protein mutation data (702 stable mutations and 2509 unstable mutations). The test set S sym contained 15 proteins, a total of 684 protein mutation data (342 stable mutations and 342 unstable mutations); the P53 dataset contained 1 protein (p53 protein), a total of 42 protein mutation data (11 stable mutations and 31 unstable mutations); the Myoglobin dataset contained 1 protein (myoglobin protein), a total of 134 protein mutation data (38 stable mutations and 96 unstable mutations). The identity of the proteins in the training set and the test set in the Protein Data Bank (Protein Data Bank), i.e. PDB ID, is shown in Table 1.

[0035] Table 1 Proteins in the training set and the test set

[0036]

[0037] The following describes the four steps of the protein mutation stability change prediction method.

[0038] Step 1: Data enhancement processing based on antisymmetry

[0039] The Q3211 dataset is an imbalanced dataset, including 702 stable mutations (ΔΔG > 0) and 2509 unstable mutations (ΔΔG < 0). Therefore, we generated a new dataset (named Q anti 3211) by data augmentation using antisymmetry. Q anti 3211 includes 2509 stable mutations and 702 unstable mutations. Finally, the balanced training set Q6422 is constructed by combining Q anti 3211 and Q3211, including 3211 stable mutations and 3211 unstable mutations.

[0040] Step 2: Extracting protein features

[0041] The sequence and structural features of each sample are extracted from the protein, including physicochemical features, secondary structure features, depth and prominence indices, evolutionary information (PSSM and SNB-PSSM), and topological features based on persistent homology. The features of each sample consider wild type, mutant (in this invention, the original protein is defined as wild type, and the mutated protein is defined as mutant), and the difference between them, which are described as follows:

[0042] (1) Physicochemical features

[0043] From the AAindex database constructed by Kawashima et al. (https: / / www.genome.jp / aaindex / ), 9 physical and chemical properties of amino acids are obtained, including the number of atoms (Atom Number, AN), hydrophobicity (Hdrpo), hydrophilicity (Hdrpi), interface propensity (Prop), isoelectric point (Isoep), mass (Mass), volume (Volume), accessible surface area (ASA), and electrostatic property classification (EPclass, positive, negative, and non-charged are represented as 1, -1, and 0, respectively), as shown in Table 2. In this way, a total of 27-dimensional features are constructed.

[0044] Table 2 9 physical and chemical properties of 20 amino acids

[0045]

[0046] (2) Secondary structure features

[0047] The secondary structure features of proteins were calculated by SPOT-1D-LM software, including three states of protein secondary structure (helix (H), strand (E) and coil (C)), eight states of protein secondary structure (3 10 helix (G), alpha-helix (H), pi-helix (I), beta-strand (E), bridge (B), turn (T), bend (S) and other (C)), relative solvent accessibility (RSA), protein backbone angles (psi angle, theta angle and tau angle), half-sphere exposure (HSE) and contact number (CN), so finally 57-dimensional features were extracted.

[0048] (3) Depth index (DPX) and protrusion index (CX)

[0049] The geometric complementarity between the mutated residue and its surrounding residues is also important for protein stability. Depth index (DPX) and protrusion index (CX) were used to represent the buried and exposed degree of atoms in protein structure, respectively. We used PSAIA (https: / / sourceforge.net / projects / psaia / ) software to calculate DPX and CX of the mutated residues, including the mean and standard deviation of DPX and CX of all atoms in the residue, as well as the mean and standard deviation of DPX and CX of the side chain atoms, finally producing 24-dimensional features.

[0050] (4) Position-specific scoring matrix (PSSM) and spatially neighboring position-specific scoring matrix

[0051] (SNB-PSSM)

[0052] PSSM gives the evolutionary conservation information of each position on the protein sequence, which can be generated by iterating PSI-BLAST (https: / / blast.ncbi.nlm.nih.gov / Blast.cgi) three times. At the same time, considering the spatial cooperativity of protein residue evolution, we previously proposed a spatially neighboring PSSM (SNB-PSSM) to extract the evolutionary information of residues. SNB-PSSM uses a spatially neighboring smoothing processing and window method to encode evolutionary information. In this calculation method, the evolutionary score of the target residue is the average of the evolutionary scores of all residues in the standard PSSM whose Cα atom distance from the target residue is not more than Figure 2 ​As shown, we extracted evolutionary scores for mutation sites using PSSM and SNB-PSSM, ultimately generating a total of 6-dimensional features.

[0053] (5) Topological features based on continuous homology

[0054] Topological features are calculated from topological invariants, which consist of independent components, loops, and cavities, and can be represented by Betti-0, Betti-1, and Betti-2, respectively. The calculation process of topological features mainly consists of two steps: the construction of simple complexes and the generation of topological features.

[0055] 5.1 Construction of Simple Complexes

[0056] In the topological space of a protein, we are concerned with the interactions between atoms. Therefore, in calculations, we will ensure that the distance to the mutant residue (Cα atom) does not exceed [a certain value]. The atoms are used to construct a point cloud (containing only three heavy atoms: C, N, and O) to extract the topological features of the mutation site. A sphere is generated centered on each point. As the sphere radius increases, when any two spheres come into contact, their centers are connected by an edge. This process is called filtering in continuous homology, and the sphere radius is the filtering radius. This process generates a series of simplicial complexes, a concept in topology referring to topological objects obtained by "gluing together" simplexes such as points, line segments, and triangles. See the appendix for details. Figure 3 .

[0057] 5.2 Generation of Topological Features

[0058] The traditional output of continuous homology is a barcode. Each individual component, loop, and cavity is represented by a bar. The start, end, and length of the bar correspond to the birth, death, and "lifetime" of the topological invariant; longer bars represent more stable topological structures. These barcode diagrams reflect the variation of topological invariants at different scales. To apply these barcodes as input to machine learning models, they need to be converted into vectors. In the barcode diagram, the horizontal axis represents the size of the filter radius, with a maximum value set to [value missing]. The entire filtration range is divided into 12 bins using the binning method. The statistics of each Bin can then be calculated. For Betti-0 bar plots, we calculate the number of bars and the number of bar death events in each Bin, resulting in 648-dimensional features. For Betti-1 and Betti-2 bar plots, we calculate the maximum, minimum, mean, standard deviation, sum of bar lengths, and the earliest birth time and the latest death time of bars, resulting in 378-dimensional features. The final integration results in a total of 1026 features. To better understand the above process, we take the protein with PDB ID 1AKK as an example to show the process of generating topological features from barcode plots (see FIG. 6). Figure 3 ).

[0059] Step 3: Training model using regression stratification sampling strategy

[0060] We designed a training scheme called regression stratification sampling to solve the problem that it is difficult to predict mutations with extreme ΔΔG values. The method consists of three steps: first, divide the range of ΔΔG values in the dataset into eight different interval intervals (labeled 1-8), including (-∞, -7], (-7, -3], (-3, -1], (-1, 0], (0, 1], (1, 3], (3, 7], (7, +∞), according to the interval of each sample, assign the corresponding label to each sample; second, divide the data in each interval into 15 equal parts, one part from each interval constitutes the validation set, and the remaining 14 parts are used as the training set; third, repeat the second step 15 times, each time using a different validation set, resulting in 15 different combinations of validation sets and training sets (see FIG. 7), then use these data to train the CNN model. Figure 4 ).

[0061] Step 4: Constructing an ensemble model based on multiple CNN sub-models

[0062] We use the 15 different sets of data obtained in step 3 for training, and choose CNN as our model framework. The network structure of CNN contains six one-dimensional convolutional layers, two dropout layers, two max pooling layers, a flatten layer and a fully connected layer. Each convolutional layer uses a rectified linear unit (ReLU) as its activation function, and the pooling size of the max pooling layer is 3. The dropout rates of the first and second dropout layers are set to 0.2 and 0.5, respectively. The batch size and the number of iterations of the CNN model are set to 20 and 1000, respectively. The Adam optimizer is used to train the CNN model with a learning rate of 0.0001. During the training process, the mean square error (MSE) is used as the loss function. Through training, 15 CNN sub-models are finally obtained, and then these 15 sub-models are integrated to build the final protein mutation stability change prediction model PMSPcnn, see Fig. 2. Figure 5 .

[0063] The training set of the PMSPcnn model is Q6422, and the three independent test sets S sym , P53, Myoglobin are used for testing. The evaluation indicators of the model are the Pearson correlation coefficient (denoted by r) and the root mean square error (denoted by σ, unit: kilocalorie / mole), and the evaluation indicator formula is as follows:

[0064]

[0065]

[0066] where x i and y i represent the sample experimental value and the predicted value (x and are their average values), respectively; n represents the total number of samples in the corresponding data set.

[0067] This invention proposes a protein mutation stability change prediction model, PMSPcnn, based on an ensemble CNN model and regression hierarchical sampling, which considers antisymmetry, spatial nearest neighbor evolutionary information, and topological features. This model can quickly and accurately predict protein stability changes caused by mutations. The model achieves higher prediction accuracy than other existing methods. This is attributed to the fact that this invention considers not only the physicochemical and structural features of residues but also residue co-evolutionary information and topological features based on persistent isohomology. Furthermore, it proposes and adopts a regression hierarchical sampling strategy. All of these factors effectively improve the prediction performance of protein mutation stability changes and can be applied to protein engineering and related drug design research. Attached Figure Description

[0068] Figure 1 Flowchart of a protein mutation stability change prediction method based on ensemble CNN model and hierarchical regression training, taking into account antisymmetry, spatial nearest neighbor evolutionary information and topological features;

[0069] Figure 2 The process of encoding the evolutionary information of the target residue in the SNB-PSSM method; (A) the spatial nearest neighbor residues (3, 4, 5, 13 and 14) of target residue 6, whose Cα atoms are located in the Cα atom of residue 6. Within the range. (B) The evolutionary information of residue 5 is defined as the average of the evolutionary information (20-dimensional vector) of its spatial nearest neighbors residues 3, 4, 5, 13 and 14 in the standard PSSM.

[0070] Figure 3 Simplex diagrams and barcodes generated for the continuous homology filtering process; (A) Simplex generated by continuous homology. (B) Continuous homology filtering process. (C) Barcode generation. Residues at mutation sites are represented by cartoons, and spheres are the residues (Cα atoms) surrounding the spheres. The area within the specified range.

[0071] Figure 4 This is a schematic diagram of the regression stratified sampling principle proposed in this invention;

[0072] Figure 5 A framework diagram of the PMSPcnn integration model;

[0073] Figure 6 The results of the PMSPcnn model in Example 1 are shown;

[0074] Figure 7 For dataset S sym Distribution of PMSCNN and ThermoNet predictions and experimental values ​​on P53 and Myoglobin. Detailed Implementation

[0075] The application will be further described in connection with the following examples, but the application is not limited to the following examples.

[0076] Example 1

[0077] The following describes the implementation process of the PMSPcnn method for predicting the mutation stability change (ΔΔG) of the protein with PDB ID 1LNI (the chain identifier of the mutation residue is A, the mutation residue number is 39, the pre-mutation residue is asparagine, and the post-mutation residue is aspartic acid). The prediction process needs to be completed under the Linux system.

[0078] 1. Preparation

[0079] Before running the protein mutation stability change prediction program, it is necessary to ensure that python 3.8 and the corresponding toolkits numpy, sklearn, xgboost, scipy and keras have been installed. Here, the PDB file of the protein 1LNI is recorded as 1LNI.pdb. The structure of the protein after mutation is obtained using Scap, and the generated PDB file is recorded as 1LNI_m.pdb. MATLAB and R software are installed, and they are used to calculate the topological features based on persistent homology. In addition, the protein mutation information is stored in an excel file, named mut_inf.csv, and the file content includes 1LNI (PDB ID), A (chain), 39 (mutation residue number), N (pre-mutation residue name) and D (post-mutation residue name).

[0080] 2. Feature extraction

[0081] (1) Physicochemical characteristics

[0082] Run the “run_phy_cha.sh” script file (command line:. / run_phy_cha.sh mut_inf.csv) to obtain the output file “PC_direct.npy”.

[0083] (2) Secondary structure characteristics

[0084] a. Download and install the “SPOT-1D-LM” program in “ https: / / github.com / jas-preet / SPOT-1D-LM.git ”;

[0085] b. Run the “run_secondary_structural.sh” script file (command line:. / run_secondary_structural.sh) to finally obtain the output file “Spot_direct.npy”.

[0086] (3) Depth index and protrusion index

[0087] a. Download and install the "psaia.exe" program from the website "https: / / sourceforge.net / projects / psaia / ";

[0088] b. Run "psaia.exe";

[0089] c. Select the “Structure Analyser” tab control, find “Analysis Types” and check “Analysis Bound” and “Analysis by Chain”, and set all parameters to default;

[0090] d. Input the PDB files “1LNI.pdb” and “1LNI_m.pdb” into the program, and then click “Run” to obtain the result files “1LNI.tbl” and “1LNI_m.tbl”;

[0091] e. Run the script file "run_dpx_cx.sh" (command line: . / run_dpx_cx.sh mut_inf.csv), and finally get the output file "psaia_direct.npy".

[0092] (4) Location-specific scoring matrix and spatial nearest neighbor location-specific scoring matrix

[0093] a. Access the server "https: / / blast.ncbi.nlm.nih.gov / Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome", enter the protein sequence of 1LNI, select "PSI-BLAST" in the "Program Selection" option, set Query Cover > 75%, E-value = 0.0001, and use the default parameters for other options. Click "multiple alignment", and after the results appear on the interface, click "download" and select 'Fasta plus gaps' to download the .fa format file and rename it to "1LNI.fa";

[0094] b. Run the script file "run_pssm.sh" (command line: . / run_pssm.sh mut_inf.csv), and finally get the output files "PSSM_direct.npy" and "PSSM_cut_direct.npy".

[0095] (5) Topological features based on continuous homology

[0096] Run the "run_topology.sh.sh" script file (command line:. / run_topology.sh mut_inf.csv) to obtain the output files "X_h0.np" and "X_h12.npy".

[0097] 3. Feature integration

[0098] All the above-mentioned feature files include "PC_direct.npy", "Spot_direct.npy", "psaia_direct.npy", "PSSM_direct.npy", "PSSM_cut_direct.npy", "X_h0.np" and "X_h12.npy". Put all the above-mentioned feature files into the folder where "feature_integrate.py" is located, and run feature_integrate.py (command line: python feature_integrate.py) to obtain the model input file "X_test_nor.npy" of the final sample.

[0099] 4. Prediction of protein stability change before and after mutation

[0100] Put "X_test_nor.npy" and "predict.py" in the same path, and run predict.py (command line: python predict.py) to obtain the final prediction result output file "Predcition_value.txt". The prediction result is printed in the bash command interface, as shown in the following figure: Figure 6 .

[0101] Example 2

[0102] The PMSPcnn model constructed in the application is tested on three independent test sets (S sym , P53, Myoglobin). The calculation process of each sample in the test set is the same as that in Example 1, and only the results are described here. On the widely used independent test set S sym , PMSPcnn is compared with 18 protein mutation stability change prediction methods developed at present, and the results are shown in Table 3.

[0103] Table 3 Comparison of performance of PMSPcnn and other methods on test set S sym

[0104]

[0105]

[0106] ​From Table 3, it is clear that PMSPcnn performs the best overall, especially on the reverse mutations, which it outperforms all other methods. PMSPcnn performs equally well on both direct and reverse mutations, with a Pearson correlation coefficient of 0.73 and a root mean square error of 1.06 kcal / mol, which reflects the unbiased nature of our model.

[0107] To verify the generalization ability of the model, we also compared our model with other models on the P53 and Myoglobin datasets. Since ThermoNet is the second best model in Table 3, we only compared our model with ThermoNet here, and the results are shown in Table 4.

[0108] Table 4. Comparison of PMSPcnn and ThermoNet on the test sets P53 and Myoglobin

[0109]

[0110] According to the results in Table 4, for the test set P53, PMSPcnn performs equally well on both direct and reverse mutations, with r and σ of 0.67 and 1.39 kcal / mol, respectively. Compared with ThermoNet, PMSPcnn improves the Pearson correlation coefficient by 49% and 20% on direct and reverse mutations, respectively. For the test set Myoglobin, PMSPcnn performs better on both direct and reverse mutations, with r and σ of 0.67 and 0.87 kcal / mol, respectively, and r and σ of 0.58 and 0.99 kcal / mol, respectively, which is much better than ThermoNet. PMSPcnn improves the Pearson correlation coefficient by 76% and 57% on direct and reverse mutations, respectively.

[0111] Figure 1 Figure 7 Figure 1 shows the distribution of predicted values by PMSPcnn and ThermoNet, as well as experimental values, on the three datasets (S sym sym ​For the dataset of PDB, the interval of experimental values is [-7.50, 7.50], the interval of predicted values of PMSPcnn is [-7.03, 6.28], and the interval of predicted values of ThermoNet is [-3.10, 2.85]; for the dataset of P53, the interval of experimental values is [-4.78, 4.78], the interval of predicted values of PMSPcnn is [-3.44, 3.42], and the interval of predicted values of ThermoNet is [-2.20, 1.88]; for the dataset of Myoglobin, the interval of experimental values is [-5.59, 5.59], the interval of predicted values of PMSPcnn is [-3.47, 3.39], and the interval of predicted values of ThermoNet is [-2.75, 2.88]. It can be seen that on the three datasets, the prediction interval of our method is the closest to the experimental value interval, which proves that the prediction ability of PMSPcnn for mutations with extreme ΔΔG values is significantly better than that of ThermoNet model, which is mainly due to the regression stratified sampling strategy we proposed.

[0112] The above results show that the newly added features (including depth and prominence index, spatial neighbor evolutionary information and topological features based on persistent homology) in the PMSPcnn method, the regression stratified sampling strategy we proposed and the ensemble model construction strategy based on multiple CNN sub-models can indeed improve the prediction performance of the model for the stability change of protein mutations, and their introduction has a good contribution to predicting the stability change of protein mutations.

Claims

1. A method for predicting the change in stability of a protein mutation based on an integrated convolutional neural network model and regression stratified sampling training, characterized in that, Considering the anti-symmetry and utilizing the spatial neighborhood evolutionary information, the persistent homology is used to extract the topological features, and the regression stratification strategy is used to train the model, which contains four steps: first, based on the anti-symmetry, the data is enhanced to balance the stable and unstable mutation sample data in the training set; second, the protein features are extracted; third, the regression stratification sampling strategy is used to train the model; and fourth, an integrated model based on multiple CNN sub-models is constructed for protein mutation stability change prediction. Before building the model, the training set and test set are first constructed: Collect experimental data on the change in stability of the protein mutation, i.e., ΔΔG: Q3421 data set; S sym data set; P53 data set; Myoglobin dataset; Q3421 as the training set, and the other three datasets S sym , P53 and Myoglobin as the test set; delete the same samples in the training set and the test set, and then use the BLAST method to remove proteins with an E-value less than 0.001 between the training set and the test set protein sequences; after processing, the final training set Q3211 contains 147 proteins, a total of 3211 protein mutation data, including 702 stable mutations and 2509 unstable mutations; the test set S sym contains 15 proteins, a total of 684 protein mutation data, including 342 stable mutations and 342 unstable mutations; the P53 dataset contains 1 p53 protein, a total of 42 protein mutation data, including 11 stable mutations and 31 unstable mutations; the Myoglobin dataset contains 1 myoglobin protein, a total of 134 protein mutation data, including 38 stable mutations and 96 unstable mutations; Step 1: Data enhancement based on anti-symmetry The Q3211 dataset is an imbalanced dataset, including 702 stable mutations and 2509 unstable mutations; a new dataset named Q anti 3211 is generated by using the antisymmetry to perform data augmentation on the data anti 3211 includes 2509 stable mutations and 702 unstable mutations; finally, the balanced training set Q6422 is constructed by combining the Q anti 3211 and Q3211 datasets, including 3211 stable mutations and 3211 unstable mutations; Step 2: Extract protein features The sequence and structure features of each sample are extracted from the protein, including physicochemical features, secondary structure features, depth and protrusion indices, evolutionary information PSSM and SNB-PSSM, and topological features based on persistent homology. The features of each sample consider the wild type, mutant type, and the difference between them, which are described as follows: (1) Physicochemical features Nine physical and chemical properties of amino acids are obtained from the AAindex database, including the number of atoms, hydrophobicity, hydrophilicity, interfacial propensity, isoelectric point, mass, volume, accessible surface area, and electrostatic properties of amino acids, a total of 27-dimensional features are constructed; (2) Secondary structure features The SPOT-1D-LM software was used to calculate the secondary structure characteristics of the protein, including three states of the secondary structure of the protein: helix (H), strand (E) and loop (C), eight states of the secondary structure of the protein: 3 10 helix (G), alpha-helix (H), pi-helix (I), beta-strand (E), bridge (B), turn (T), bend (S) and other (C), relative solvent accessibility, protein backbone angles: psi angle, phi angle, theta angle and tau angle, hemisphere exposure and contact number, finally 57-dimensional features were extracted; (3) Depth index DPX and protrusion index CX Depth index and protrusion index are used to represent the buried and exposed degree of atoms in protein structure, PSAIA software is used to calculate DPX and CX of mutant residues, including the mean and standard deviation of DPX and CX of all atoms in the residue, as well as the mean and standard deviation of DPX and CX of side chain atoms, finally 24-dimensional features are generated; (4) Position-specific scoring matrix PSSM and spatial neighborhood position-specific scoring matrix SNB-PSSM PSSM gives the evolutionary conservation information of each position in the protein sequence, which can be generated by iterating PSI-BLAST three times, and considering the spatial coordination of protein residue evolution, SNB-PSSM is used to extract the evolutionary information of residues based on spatial neighborhood, SNB-PSSM uses smoothing processing and window method based on spatial neighborhood to encode evolutionary information, the evolutionary score of the target residue is the average of the evolutionary scores of all residues within a distance of 7.5 Å from the Cα atom of the target residue, only the evolutionary scores of the mutant site are extracted by PSSM and SNB-PSSM, finally 6-dimensional features are generated; (5) Topological features based on persistent homology Topological features are calculated from topological invariants, which are composed of independent groups, loops and cavities, represented by Betti-0, Betti-1 and Betti-2 respectively; The calculation process of topological features consists of two steps: construction of simplicial complex and generation of topological features; 5.1 Construction of simplicial complex In the protein topology space, atoms within a distance of 12 Å from the C a atom of the mutated residue are constructed into a point cloud containing only C, N and O atoms, a sphere is generated with each point in the point cloud as the center, the radius of the sphere is the filtering radius, a series of simple complexes are generated, and the simple complexes are topological objects obtained by bonding points, line segments and triangular simplices; 5.2 Generation of topological features The traditional output of persistent homology is a bar code, and each independent component, loop and cavity is represented by a bar Bar. In the bar code, the horizontal axis represents the size of the filtering radius, and the maximum value is set to 12 Å. The entire filtering range is divided into 12 bins of 1.0 Å using the binning method, and then the statistics of each bin are calculated. For the Betti-0 bar code, the number of Bar bars and the number of Bar death events in each bin are calculated, and a total of 648 features are obtained. For the Betti-1 and Betti-2 bar codes, the maximum, minimum, mean, standard deviation, and total of the Bar length are calculated, as well as the initial birth time and the final death time of the Bar, and a total of 378 features are obtained. Finally, a total of 1026 features are integrated. Step 3: Train the model using a regression stratified sampling strategy A regression stratified sampling training scheme is designed, which consists of three steps: first, divide the range of ΔΔG values in the data set into eight different interval intervals, including (-∞, -7], (-7, -3], (-3, -1], (-1, 0], (0, 1], (1, 3], (3, 7], (7, +∞), according to the interval of each sample, assign the corresponding label to each sample; second, divide the data in each interval into 15 equal parts, one part from each interval forms a validation set, and the remaining 14 parts are used as a training set; third, repeat step two 15 times, each time using a different validation set, to obtain 15 different combinations of validation sets and training sets, and then train the CNN model; Step 4: Build an integrated model based on multiple CNN sub-models Train using the 15 different sets of data obtained in step 3, and select CNN as the model framework; the network structure of CNN includes six one-dimensional convolutional layers, two random deactivation Dropout layers, two maximum pooling MaxPooling1D layers, a flattening Flatten layer and a fully connected Fully Connected layer, each convolutional layer uses a linear rectifier function ReLU as its activation function, and the pooling size of the maximum pooling layer is 3; the dropout rates of the first and second Dropout layers are set to 0.2 and 0.5 respectively, the batch size and iteration number of the CNN model are set to 20 and 1000 respectively, and the Adam optimizer is used to train the CNN model with a learning rate of 0.0001, during the training process, the mean square error is used as the loss function, and finally 15 CNN sub-models are obtained by training, then the 15 sub-models are integrated to build the final protein mutation stability change prediction model PMSPcnn.

2. The method according to claim 1, characterized in that We construct a balanced training set Q6422 by augmenting the data using anti-symmetry and test the model on three independent test sets S sym , P53, Myoglobin, and the evaluation metrics are the Pearson correlation coefficient r and the root mean square error Unit: kcal / mol, and the evaluation metrics are as follows: （1） （2） where x i and y i represent the experimental and predicted values, respectively, is the average of x i , is the average of y i ; and n represents the total number of samples in the corresponding data set.

Images