Protein mutation stability change prediction method based on integrated convolutional neural network model and regression hierarchical training

Abstract

The invention discloses a protein mutation stability change prediction method based on an integrated convolutional neural network model and regression hierarchical training, and belongs to the technical field of protein stability change prediction. Firstly, a training set and a test set are constructed, and four data sets are collected and arranged. The method comprises the following four steps: 1, performing enhancement processing on data based on antisymmetry, and balancing and training stable and unstable abrupt sample data in a set; 2, extracting protein characteristics; 3, training the model by adopting a regression stratified sampling strategy; and 4, constructing an integrated model based on a plurality of CNN sub-models for predicting the stability change of the protein mutation. Spatial neighbor evolution information is introduced for the first time; a regression stratified sampling strategy is proposed for the first time and is successfully applied to the training process of the model, so that the prediction capability of the model for the extreme delta delta G is effectively improved; the trained CNN sub-models are combined, information of all samples in a training set is fully utilized, and the generalization ability of the model is enhanced.

Description

Technical Field

[0001] The invention 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 Art

[0002] Many human genetic diseases are directly associated with missense mutations in proteins [1] Missense mutations often lead to an increase (stabilizing mutation) or decrease (destabilizing mutation) in protein stability, thus affecting its function and causing disease. [2] Currently, accurately predicting the stability changes caused by protein mutations is an important and challenging task. The protein stability changes caused by mutations are generally quantified by measuring the change in folding free energy (ΔΔG) before and after protein mutation. [1] Currently, experimental methods such as differential scanning calorimetry and thermofluorescence spectroscopy can detect changes in protein stability. [3] However, these methods are time-consuming and cumbersome. Therefore, there is an urgent need to develop an efficient computational method to predict protein stability changes caused by mutations.

[0003] Several theoretical methods have been developed for predicting changes in protein stability (ΔΔG), which can be divided into three main categories: energy-based, knowledge-based, and machine learning-based methods. [4] Energy-based methods include free energy perturbation (FEP) [5] , which is an early developed method based on molecular dynamics simulation and is generally considered to have high accuracy. The FoldX model is a traditional energy-based method that uses empirical force fields and energy terms to calculate the effects of mutations 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 environment-specific substitution tables (ESSTs) with conformational constraints to calculate changes in protein stability before and after mutation, while DDGun only used evolutionary information to calculate differences in stability. Both tools use statistical methods with limited generalization capabilities. 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] Etc. These methods use traditional tree models or deep neural networks and take into account various features, including sequence, structure, and physicochemical characteristics, which improves the generalization ability of the model to a certain extent.

[0004] Most methods for predicting changes in protein mutation stability generally have two problems: first, the prediction ability for stable mutations is poor; second, the prediction ability for mutations with extreme ΔΔG values ​​(highly unstable or highly stable mutations) is insufficient. 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 dataset are unstable).

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

[16] . We consider using this property to construct a balanced training set to improve the prediction bias caused by unbalanced data. For the second problem, the main reason is that there are generally few mutation samples 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 small samples during the training process is important for the model to accurately predict extreme values. As far as we know, existing methods have ignored this point, so there are currently no measures to solve this problem. It is worth noting that mutations with extreme ΔΔG values ​​are generally more likely to have a greater impact on the structure of the protein and thus lead to the occurrence of diseases. In order to solve this problem, we designed a regression stratified sampling training method, which will help improve the model's predictive ability for extreme ΔΔG values.

[0005] In addition, feature extraction is also crucial for the model. We considered physicochemical features, Position Specific Scoring Matrix (PSSM), secondary structure features, as well as depth index and prominence index. In addition, we recently proposed a new method for encoding evolutionary information, namely Spacial Neighbor-Based Position-Specific Scoring Matrix (SNB-PSSM)

[17] , which takes into account the characteristics of residue space co-evolution (averaging neighboring residues in space) when encoding residue evolution information, and has achieved good results in protein-RNA binding site prediction. In addition to the above traditional features, we believe that further consideration of topological structure 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 in many aspects, such as protein-ligand binding affinity prediction, protein interaction prediction and drug virtual screening [19-21] Therefore, we believe that topological features based on persistent coherence will also help improve the model's ability to predict changes in protein stability.

[0006] In recent years, deep learning has been widely used in the field of protein structure function prediction, such as protein structure prediction, functional annotation, interaction prediction and drug design. Among them, convolutional neural networks (CNN) are more efficient and accurate in processing large-scale protein data, and have advantages such as translation invariance and parameter sharing, and perform well in capturing local features in the data. Most deep learning will look for the model with the best performance as the final prediction model during the training process, but the samples used as the validation set in the training do not participate in the training of the model, resulting in the loss of information contained in this part of the samples. Ensemble Method

[22] This problem can be solved well. Therefore, we consider using an integrated method to build a CNN model framework.

[0007] In summary, we consider antisymmetry to solve the problem of data imbalance; we extract protein sequence and structural features, and introduce spatial neighbor evolution information and topological features based on continuous homology to fully consider the features related to protein stability; we propose a stratified sampling strategy, called regression stratification cross-validation (RScv), to solve the problem of poor prediction of mutations with extreme ΔΔG values; and we propose an ensemble method to combine the trained CNN sub-models to make full use of the information of all samples in the training set. Finally, we constructed an unbiased model of protein mutation stability changes based on an integrated 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 Mutationson 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: toward 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 ofcomputational tools to predict the effects of mutations on protein stabilityin the absence of experimental structures[J].Brief Bioinform,2022,23(2):bbac025.

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

[0014] [6]GUEROIS R,NIELSEN J E,SERRANO L.Predicting changes in thestability of proteins and protein complexes:a study of more than1000mutations[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 forpredicting 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 methodfor the prediction of protein stability changes upon single and multiplepoint variations[J].BMC Bioinformatics,2019,20(Suppl 14):335.

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

[0018]

[10] PIRESD E,ASCHER D B,BLUNDELL T L.mCSM:predicting the effects ofmutations 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 predictingeffects of mutations on protein stability using an integrated computationalapproach[J].Nucleic Acids Res,2014,42(Web Server issue):W314-W319.

[0020]

[12] CAPRIOTTI E,FARISELLIP,CASADIO R.I-Mutant2.0:predicting stabilitychanges upon mutation from the protein sequence or structure[J].Nucleic AcidsRes,2005,33(Web Server issue):W306-W310.

[0021]

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

[0022]

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

[0023]

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

[0024]

[16] THILTGEN G,GOLDSTEIN R A.Assessing predictors of changes inprotein 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 PSSMused 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 of the Invention

[0031] The purpose of the present invention is to predict the mutation stability change (ΔΔG) when an amino acid mutation occurs at a certain site in a known protein, which is helpful for the study of protein structure-function relationship and can provide important information for protein engineering and drug design. The characteristics of this method are: ① Based on antisymmetry, the data is enhanced and processed to finally generate a balanced data set (containing the same number of stable and unstable mutations), which is conducive to solving the problem of prediction bias for unstable mutation samples and improving the ability to predict stable mutations. ② In terms of feature extraction, in addition to traditional protein sequence and structural features, the present invention also introduces a new evolutionary information encoding method for the first time, namely, a position-specific scoring matrix based on spatial neighbors, which was previously developed by our group (J Mol Recognit. 2021Jun 1; 6(34): e2887). In addition, topological features based on continuous homology are also used to predict protein stability changes. These features can play an important role in predicting protein stability. ③ For the first time, a regression stratified sampling strategy is proposed and applied to the prediction task of protein mutation stability changes, which effectively improves the model's prediction ability for extreme ΔΔG. ④Different from previous methods that use deep learning to predict changes in protein mutation stability, we proposed a new and effective integrated model based on an integrated approach that contains multiple CNN sub-models, which is conducive to enhancing the generalization ability of the model.

[0032] The present invention is based on an ensemble convolutional neural network model and a regression stratified sampling strategy, and is a protein mutation stability change prediction method (anensemble convolutional neural network model with regression stratification cross-validation for Protein Stability changes upon Mutations Prediction, PMSPcnn) that considers antisymmetry and utilizes spatial neighbor evolution information and continuous homology to extract topological features. The method comprises four steps: first, enhancing the data based on antisymmetry to balance the sample data of stable and unstable mutations in the training set; second, extracting protein features; third, adopting a regression stratified sampling strategy to train the model; and fourth, constructing an integrated model based on multiple CNN sub-models for protein mutation stability change prediction (such as the attached Figure 1 ).

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

[0034] Collect experimental data on protein mutation stability changes (ΔΔG) from recently published literature: Q3421 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 the training set, and the other three datasets S sym , P53 and Myoglobin are the test set. The samples that are the same in the training set and the test set are deleted, and then the BLAST method (https: / / ftp.ncbi.nlm.nih.gov / blast / executables / blast+ / LATEST / ) is used to remove proteins whose E-value (an indicator for evaluating the similarity of two sequences, the smaller the E-value, the higher the similarity of the two sequences) of the protein sequence in the training set and the test set is less than 0.001. After processing, the training set (named Q3211) contains 147 proteins and a total of 3211 protein mutation data (702 stable mutations and 2509 unstable mutations). Test set S sym The dataset contains 15 proteins, with a total of 684 protein mutation data (342 stable mutations and 342 unstable mutations); the P53 dataset contains 1 protein (p53 protein), with a total of 42 protein mutation data (11 stable mutations and 31 unstable mutations); the Myoglobin dataset contains 1 protein (myoglobin protein), with a total of 134 protein mutation data (38 stable mutations and 96 unstable mutations). The identification numbers (Identity) of the proteins in the training set and the test set in the Protein Data Bank (PDB ID) are shown in Table 1.

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

[0036]

[0037] The four steps of the method for predicting changes in protein mutation stability are described below.

[0038] Step 1: Enhance the data based on antisymmetry

[0039] The Q3211 dataset is an unbalanced dataset, including 702 stable mutations (ΔΔG>0) and 2509 unstable mutations (ΔΔG<0). Therefore, we used antisymmetry to perform data augmentation on the data and generated a new dataset (named Q anti 3211).Q anti 3211 includes 2509 stable mutations and 702 unstable mutations. Finally, the constructed balanced training set Q6422 is composed of Q anti The 3211 and Q3211 datasets are combined, including 3211 stable mutations and 3211 unstable mutations.

[0040] Step 2: Extract protein features

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

[0042] (1) Physical and chemical characteristics

[0043] Nine physicochemical properties of amino acids were obtained from the AAindex database constructed by Kawashima et al. (https: / / www.genome.jp / aaindex / ), including the number of atoms (AN), hydrophobicity (Hdrpo), hydrophilicity (Hdrpi), interface tendency (Propensity), isoelectric point (Isoep), mass, volume, accessible surface area (ASA), and electrostatic property classification (EPclass, positive, negative, and non-charged are represented by 1, -1, and 0, respectively) of each amino acid, see Table 2. In this way, a total of 27 dimensions of features were constructed.

[0044] Table 2 9 physicochemical property parameters of 20 amino acids

[0045]

[0046] (2) Secondary structure characteristics

[0047] SPOT-1D-LM software was used to calculate the secondary structural features of proteins, including three states of protein secondary structure (helix (H), chain (E) and loop (C)), eight states of protein secondary structure (3 10 -helix (G), α-helix (H), π-helix (I), β-strand (E), bridge (B), turn (T), bend (S) and other (C)), relative solvent accessibility (Relative Solvent Accessibility, RSA), protein main chain angle (ψ angle, angle, θ angle and τ angle), hemisphere exposure (Half-Sphere Exposure, HSE) and contact number (Contact Number, CN), so that a total of 57 dimensions of features are extracted.

[0048] (3) Depth Index (DPX) and Prominence Index (CX)

[0049] The geometric complementarity between the mutated residues of a protein and the surrounding residues is also important for protein stability. The depth index (DPX) and the protrusion index (CX) are used to characterize the degree of burial and exposure of atoms in the protein structure, respectively. We use the PSAIA (https: / / sourceforge.net / projects / psaia / ) software to calculate the 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 side chain atoms, and finally generate a total of 24-dimensional features.

[0050] (4) Position-specific scoring matrix (PSSM) and spatial neighbor position-specific scoring matrix

[0051] (SNB-PSSM)

[0052] PSSM gives the evolutionary conservation information of each position in the protein sequence and can be generated by iterating PSI-BLAST (https: / / blast.ncbi.nlm.nih.gov / Blast.cgi) three times. At the same time, considering the spatial coordination of protein residue evolution, the spatial neighbor-based PSSM (SNB-PSSM) we proposed earlier is also used to extract the evolutionary information of residues. SNB-PSSM uses spatial neighbor-based smoothing and window methods to encode evolutionary information. In this calculation method, the evolutionary score of the target residue is the Cα atom distance of the target residue in the standard PSSM that is no more than The average evolutionary score of all residues in Figure 2As shown in Figure 2, we use PSSM and SNB-PSSM to extract only the evolutionary scores of mutation sites, and finally generate a total of 6-dimensional features.

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

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

[0055] 5.1 Construction of simplicial complex

[0056] In protein topological space, we care about the interaction between atoms. In this way, when calculating, we will be no more than 100% away from the mutated residue (Cα atom) The atoms are constructed into 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 with each point as the center. As the radius of the sphere increases, when any two spheres touch each other, their centers are connected by an edge. This process is called the filtering process in continuous homology, and the radius of the sphere is the filtering radius. In this process, a series of simplicial complexes will be generated. It is a concept in topology, which refers to the topological objects obtained by "gluing" simplexes such as points, line segments, and triangles. The specific process is shown in the attached Figure 3 .

[0057] 5.2 Generation of topological features

[0058] The traditional output of persistent homology is a barcode. Each independent component, ring, and cavity is represented by a bar. The start, end, and length of the bar correspond to the birth, death, and "survival time" of the topological invariant, and longer bars represent more stable topological structures. These barcodes reflect the changes in topological invariants at different scales. In order to apply these barcodes to the input of machine learning models, they need to be converted into vectors. In the bar chart, the horizontal axis represents the size of the filter radius, and the maximum value is set to The entire filtering range is divided into 12 bins using the Binning Method. Bins, and then the statistics of each Bin can be calculated. For the Betti-0 bar chart, we count the number of bars and the number of bar death events in each Bin, and get a total of 648-dimensional features. For the Betti-1 and Betti-2 bar charts, we calculate the maximum, minimum, mean, standard deviation, sum of the Bar length, as well as the Bar's initial birth time and final death time, and get a total of 378-dimensional features. Finally, a total of 1026 features were integrated. In order to better understand the above process, we take the protein with PDB ID 1AKK as an example to demonstrate the process of generating topological features from the barcode map (see Appendix Figure 3 ).

[0059] Step 3: Train the model using regression stratified sampling strategy

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

[0061] Step 4: Build an integrated model based on multiple CNN sub-models

[0062] We use 15 different sets of data obtained in step 3 for training and select CNN as our model framework. The network structure of CNN contains six one-dimensional convolutional layers, two random inactivation (Dropout) layers, two maximum pooling (MaxPooling1D) layers, a flatten (Flatten) layer and a fully connected (Fully Connected) layer. Each convolutional layer uses the linear rectified unit (ReLU) as its activation function, and the pooling size of the maximum pooling layer is 3. The loss rates of the first and second dropout layers are set to 0.2 and 0.5, respectively. The batch size (batch size) and the number of iterations (epochs) 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. After training, 15 CNN sub-models are finally obtained, and then these 15 sub-models are integrated to construct the final protein mutation stability change prediction model PMSPcnn, see the attached Figure 5 .

[0063] The PMSPcnn model training set uses Q6422 and is tested on three independent test sets S sym , P53, and Myoglobin. The evaluation indicators of the model are the Pearson correlation coefficient (expressed by r) and the root mean square error (expressed by σ, unit kcal / mol). The evaluation indicator formula is as follows:

[0064]

[0065]

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

[0067] The present invention proposes a protein mutation stability change prediction model PMSPcnn based on an integrated CNN model and regression stratified sampling that considers antisymmetry, spatial neighbor evolution information and topological features, and can quickly and accurately predict protein stability changes caused by mutations. The model has a higher accuracy rate in prediction performance than other existing prediction methods. The reason is that the present invention not only considers the physicochemical and structural characteristics of the residues, but also considers the co-evolutionary information of the residues, as well as the topological features based on continuous homology. In addition, a regression stratified sampling strategy is proposed and adopted. All of these effectively improve the prediction effect of protein mutation stability changes and can be used for protein engineering and related drug design research. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 Flowchart of the method for predicting protein mutation stability changes based on an integrated CNN model and regression hierarchical training, taking into account antisymmetry, spatial neighbor evolution information, and topological features;

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

[0070] Figure 3 Simplex graphs and barcodes generated by the continuous homology filtering process; (A) Simplex generated by continuous homology. (B) Continuous homology filtering process. (C) Generation of barcodes. The residue at the mutation site is represented by a cartoon, and the spheres are the residues (Cα atoms) around the residues. Area within the range.

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

[0072] Figure 5 This is the framework diagram of the PMSPcnn integrated model;

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

[0074] Figure 7 For the dataset S sym Distribution of PMSPcnn and ThermoNet predicted values ​​and experimental values ​​on P53 and Myoglobin. DETAILED DESCRIPTION

[0075] The present invention will be further described below in conjunction with examples, but the present invention is not limited to the following examples.

[0076] Example 1

[0077] The following is an introduction to the implementation process of the PMSPcnn method to predict the mutation stability change (ΔΔG) of a protein with PDB ID 1LNI (the chain identifier of the mutated residue is A, the mutated residue number is 39, the residue before mutation is asparagine, and the residue after mutation 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, make sure that python 3.8 and the corresponding toolkits numpy, sklearn, xgboost, scipy and keras are installed. Here, we record the PDB file corresponding to protein 1LNI as 1LNI.pdb. Use Scap to obtain the structure of the protein after mutation, and the generated PDB file is recorded as 1LNI_m.pdb. Install MATLAB and R software and use them to calculate topological features based on continuous homology. In addition, store the protein mutation information in an excel file named mut_inf.csv, and the file content includes 1LNI (PDB ID), A (chain), 39 (mutation residue number), N (residue name before mutation) and D (residue name after mutation).

[0080] 2. Feature extraction

[0081] (1) Physical and chemical characteristics

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

[0083] (2) Secondary structure characteristics

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

[0085] b. Run the "run_secondary_structural.sh" script file (command line: . / run_secondary_structural.sh), and finally obtain the output file "Spot_direct.npy".

[0086] (3) Depth Index and Prominence 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 "Analyse Bound" and "Analyse by Chain", and all parameters are set 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 "run_dpx_cx.sh" script file (command line: . / run_dpx_cx.sh mut_inf.csv), and finally obtain the output file "psaia_direct.npy".

[0092] (4) Position-specific scoring matrix and spatial neighbor position-specific scoring matrix

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

[0094] b. Run the "run_pssm.sh" script file (command line: . / run_pssm.sh mut_inf.csv), and finally obtain 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) and obtain the output files "X_h0.np" and "X_h12.npy".

[0097] 3. Feature Integration

[0098] All the feature files obtained above 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 feature files mentioned above into the folder where "feature_integrate.py" is located, run feature_integrate.py (command line: python feature_integrate.py), and get the model input file "X_test_nor.npy" of the final sample.

[0099] 4. Predict protein stability changes before and after mutation

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

[0101] Example 2

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

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

[0104]

[0105]

[0106] It can be clearly seen from Table 3 that PMSPcnn has the best overall performance among all methods, especially in predicting reverse mutations. PMSPcnn has the same performance in predicting direct mutations and reverse mutations, with a Pearson correlation coefficient of 0.73 and a root mean square error of 1.06 kcal / mol, which reflects that our model has good unbiasedness.

[0107] In order to verify the generalization ability of the model, we compared it with other models on the P53 and Myoglobin datasets. Since ThermoNet is second only to our model in Table 3, we only compare it with the ThermoNet model here. The results are shown in Table 4.

[0108] Table 4 Comparison of performance between PMSPcnn and ThermoNet on test sets P53 and Myoglobin

[0109]

[0110] According to the results in Table 4, for the test set P53, PMSPcnn has the same prediction performance for direct mutations and reverse mutations, with r and σ being 0.67 and 1.39 kcal / mol, respectively. Compared with ThermoNet, PMSPcnn's Pearson correlation coefficients for direct mutations and reverse mutations are improved by 49% and 20%, respectively; for the test set Myoglobin, PMSPcnn's r and σ are 0.67 and 0.87 kcal / mol for direct mutations, and 0.58 and 0.99 kcal / mol for reverse mutations, respectively, which are much better than ThermoNet's prediction results, and the Pearson correlation coefficients for direct mutations and reverse mutations are improved by 76% and 57%, respectively.

[0111] Attached Figure 7 The results are shown on three datasets (S sym , P53, Myoglobin) PMSPcnn and ThermoNet predicted values ​​and experimental values. symFor the P53 dataset, the experimental value interval is [-7.50, 7.50], the predicted value interval of PMSPcnn is [-7.03, 6.28], and the predicted value interval of ThermoNet is [-3.10, 2.85]. For the Myoglobin dataset, the experimental value interval is [-5.59, 5.59], the predicted value interval of PMSPcnn is [-3.47, 3.39], and the predicted value interval of ThermoNet is [-2.75, 2.88]. It can be seen that on the three data sets, our prediction intervals are the closest to the experimental value intervals, proving that PMSPcnn's prediction ability for mutations with extreme ΔΔG values ​​is significantly better than that of the ThermoNet model, which is mainly due to the regression stratified sampling strategy we proposed.

[0112] The above results show that the newly added features in the PMSPcnn method (including depth and prominence index, spatial neighbor evolution information and topological features based on continuous homology), the regression stratified sampling strategy we proposed and the integrated model construction strategy based on multiple CNN sub-models can indeed improve the model's prediction performance for changes in protein mutation stability. Their introduction has made a good contribution to the prediction of changes in protein mutation stability.

Claims

1. A method for predicting protein mutation stability changes based on an integrated convolutional neural network model and regression stratified sampling training, characterized in that: Considering antisymmetry and utilizing spatial neighbor evolution information and continuous homology to extract topological features, a regression hierarchical strategy was used to train the model, which included four steps: first, data enhancement was performed based on antisymmetry to balance the sample data of stable and unstable mutations in the training set; second, protein features were extracted; third, the model was trained using a regression hierarchical sampling strategy; and fourth, an integrated model based on multiple CNN sub-models was constructed to predict changes in protein mutation stability.

2. The method for predicting protein mutation stability changes based on an integrated convolutional neural network model and regression stratified sampling training according to claim 1, characterized in that: Before building the model, we first constructed the training set and test set: Collect experimental data on protein mutation stability changes, i.e., ΔΔG, from recently published literature: Q3421 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 the training set, and the other three datasets S sym , P53 and Myoglobin as the test set; delete the samples in the training set that are the same as those in the test set, and then use the BLAST method (https: / / ftp.ncbi.nlm.nih.gov / blast / executables / blast+ / LATEST / ) to remove the proteins whose protein sequences in the training set and the test set have an E value less than 0.001; after processing, the training set (named Q3211) contains 147 proteins and a total of 3211 protein mutation data (702 stable mutations and 2509 unstable mutations); the test set S sym Contains 15 proteins, with a total of 684 protein mutation data (342 stable mutations and 342 unstable mutations); P53 dataset contains 1 protein (p53 protein), with a total of 42 protein mutation data (11 stable mutations and 31 unstable mutations); Myoglobin dataset contains 1 protein (myoglobin protein), with a total of 134 protein mutation data (38 stable mutations and 96 unstable mutations); Step 1: Enhance the data based on antisymmetry The Q3211 dataset is an unbalanced dataset, including 702 stable mutations and 2509 unstable mutations. We used antisymmetry to perform data augmentation on the data to generate a new dataset named Q anti 3211,Q anti 3211 includes 2509 stable mutations and 702 unstable mutations; Finally, the constructed balanced training set Q6422 consists of Q anti The 3211 and Q3211 datasets were combined, including 3211 stable mutations and 3211 unstable mutations; Step 2: Extract protein features Sequence and structural features of each sample are extracted from the protein, including physicochemical features, secondary structure features, depth and protrusion index, evolutionary information (PSSM and SNB-PSSM), and topological features based on continuous homology. The features of each sample consider the wild type, mutant type (the original protein is the wild type, and the mutated protein is the mutant type), and the difference between them, as described below: (1) Physical and chemical characteristics Nine physicochemical properties of amino acids were obtained from the AAindex database constructed by Kawashima et al. (https: / / www.genome.jp / aaindex / ), including the atomic number (AN), hydrophobicity (Hdrpo), hydrophilicity (Hdrpi), interface tendency (Propensity), isoelectric point (Isoep), mass, volume, accessible surface area (ASA) and electrostatic property class (EPclass, positive, negative and non-charged are represented by 1, -1 and 0, respectively) of each amino acid, as shown in Table 2. A total of 27 dimensions of features were constructed. Table 2 9 physicochemical property parameters of 20 amino acids (2) Secondary structure characteristics SPOT-1D-LM software was used to calculate the secondary structural features of proteins, including three states of protein secondary structure (helix (H), chain (E) and loop (C)), eight states of protein secondary structure (3 10 -helix (G), α-helix (H), π-helix (I), β-strand (E), bridge (B), turn (T), bend (S) and other (C)), relative solvent accessibility (Relative Solvent Accessibility, RSA), protein main chain angle (ψ angle, angle, θ angle and τ angle), hemisphere exposure (Half-SphereExposure, HSE) and contact number (Contact Number, CN), and finally a total of 57 dimensions of features were extracted. (3) Depth Index (DPX) and Prominence Index (CX) The geometric complementarity between the mutated residues of a protein and the surrounding residues is also important for protein stability. The depth index (DPX) and the protrusion index (CX) are used to characterize the degree of burial and exposure of atoms in the protein structure, respectively. The PSAIA (https: / / sourceforge.net / projects / psaia / ) software is used to calculate the 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 side chain atoms, and finally a total of 24-dimensional features are generated; (4) Position-specific scoring matrix (PSSM) and spatial neighbor position-specific scoring matrix (SNB-PSSM) PSSM gives the evolutionary conservation information of each position in the protein sequence. It can be generated by iterating PSI-BLAST (https: / / blast.ncbi.nlm.nih.gov / Blast.cgi) three times. Taking into account the spatial coordination of protein residue evolution, PSSM based on spatial neighbors (SNB-PSSM) is used to extract the evolutionary information of residues. SNB-PSSM uses smoothing and windowing methods based on spatial neighbors to encode evolutionary information. In this calculation method, the evolutionary score of the target residue is the Cα atom distance of the target residue in the standard PSSM that is no more than 1 / 3 of the target residue. The average evolutionary scores of all residues are taken, and the evolutionary scores of only the mutation sites are extracted through PSSM and SNB-PSSM, and finally a total of 6-dimensional features are generated; (5) Topological features based on continuous homology The topological features are calculated from the topological invariants, which consist of independent components, rings, and cavities, which are denoted by Betti-0, Betti-1, and Betti-2, respectively; The calculation process of topological features is mainly divided into two steps: construction of simplicial complex and generation of topological features; 5.1 Construction of simplicial complex In protein topological space, we are concerned about the interaction between atoms. When calculating, we will consider the distance from the mutant residue (Cα atom) to be no more than The atoms are constructed into 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 with each point as the center. As the radius of the sphere increases, when any two spheres touch each other, their centers will be connected by an edge. This process is called the filtering process in continuous homology. The radius of the sphere is the filtering radius. In this process, a series of simplicial complexes will be generated. It is a concept in topology, which refers to the topological objects obtained by "gluing" simplexes such as points, line segments, and triangles. 5.2 Generation of topological features The traditional output of continuous homology is a barcode. Each independent component, ring, and cavity is represented by a bar. The start, end, and length of the bar correspond to the birth, death, and "survival time" of the topological invariant. Longer bars represent more stable topological structures. These barcodes reflect the changes in topological invariants at different scales. In order to apply these barcodes to the input of machine learning models, they need to be converted into vectors. In the bar chart, the horizontal axis represents the size of the filter radius, and the maximum value is set to The entire filtering range is divided into 12 bins using the Binning Method. Bins, and then the statistics of each Bin can be calculated; for the Betti-0 bar chart, the number of bars and the number of bar death events in each Bin are calculated, and a total of 648-dimensional features are obtained; For Betti-1 and Betti-2 bar charts, the maximum, minimum, mean, standard deviation, sum, and the initial birth time and final death time of the bar length are calculated, resulting in a total of 378 features; Finally, a total of 1026 features were obtained through integration; Step 3: Train the model using regression stratified sampling strategy A regression stratified sampling training scheme is designed. The method is divided into three steps: the first step is to divide the ΔΔG value range in the data set into eight intervals with different intervals (marked as 1-8), including (-∞, -7], (-7, -3], (-3, -1], (-1, 0], (0, 1], (1, 3], (3, 7], (7, +∞), and assign the corresponding label to each sample according to the interval in which each sample is located; the second step is to divide the data in each interval into 15 parts, one of which from each interval constitutes the validation set, and the remaining 14 parts are used as training sets; the third step is to repeat the second step 15 times, each time using a different validation set, to obtain 15 sets of validation sets and training sets with different combinations, and then use these data to train the CNN model; Step 4: Build an integrated model based on multiple CNN sub-models The 15 different sets of data obtained in step 3 were used for training, and CNN was selected as the model framework. The network structure of CNN included six one-dimensional convolutional layers, two random dropout layers, two maximum pooling (MaxPooling1D) layers, a flatten layer and a fully connected (Fully Connected) layer. Each convolutional layer used the rectified linear unit (ReLU) as its activation function, and the pooling size of the maximum pooling layer was 3. The loss rates of the first and second dropout layers were set to 0.2 and 0.5, respectively. The batch size and number of iterations (epochs) of the CNN model were set to 20 and 1000, respectively. The Adam optimizer was used to train the CNN model with a learning rate of 0.0001. During the training process, the mean square error (MSE) was used as the loss function. 15 CNN sub-models were finally obtained through training, and then these 15 sub-models were integrated to construct the final protein mutation stability change prediction model PMSPcnn.

3. The method according to claim 2, characterized in that Before building the model, we first constructed a training set and a test set. The identification numbers (Identity) of the proteins in the training set and the test set in the Protein Data Bank (PDB ID) are shown in Table 1. Table 1 Proteins in the training and test sets 4. The method according to claim 2, characterized in that We use antisymmetry to enhance the data to construct a balanced training set Q6422 and three independent test sets S sym , P53, and Myoglobin. The evaluation indicators of the model are the Pearson correlation coefficient (expressed by r) and the root mean square error (expressed by σ, unit kcal / mol). The evaluation indicator formula is as follows: where x i and i represent the sample experimental value and the predicted value (x and is their average value); n represents the total number of samples in the corresponding data set.

Images