Protein sequence design method based on biomolecular interaction domain enhancement

Abstract

The application discloses a protein sequence design method based on biomolecular interaction domain enhancement, which comprises the following steps: inputting a protein main chain skeleton three-dimensional coordinate information of a size of L*N*3 to be subjected to sequence design; obtaining a protein sequence in contact with a biomolecule and an interaction domain interval; clustering the obtained sequence and taking a representative sequence of each cluster as a training set; extracting three-dimensional structure, secondary structure, solvent accessibility and functional annotation feature representation of each training sample; using a LoRA algorithm to fine-tune the last ten transformer modules of a general multi-modal protein language model ESM3, and giving greater weight to the loss of a mask residue located in the interaction domain interval; and inputting the atomic coordinates of the protein main chain skeleton to be subjected to sequence design into the trained model to obtain a target sequence. On the one hand, the application utilizes multi-modal information of a large amount of proteins; and on the other hand, the application can generate a more robust and reasonable functional protein sequence.

Description

Technical Field

[0001] The present invention relates to the fields of structural biology, artificial intelligence and computer application technology, and more specifically, to a protein sequence design method based on enhanced biomolecular interaction domains. Background Art

[0002] Protein sequence design involves designing a primary amino acid sequence that fits the protein's backbone topology, starting from the protein's backbone. Because artificially designed proteins have human-defined functions, they have broad and profound implications for the life sciences. For example, designing specific protein sequences can develop new therapeutic drugs or enhance the efficacy of existing drugs; designing more stable and efficient enzymes can improve the efficiency and reduce costs of industrial processes; and it can be used to design highly immunogenic antigens capable of eliciting effective immune responses. Therefore, precise protein sequence design is crucial for life science fields such as drug discovery, vaccine development, and enzyme engineering. With the rise of artificial intelligence (AI), the de novo design of new protein sequences using deep learning has garnered widespread attention.

[0003] At present, the methods specifically used for protein sequence design include: ProteinMPNN (Dauparas J, Anishchenko I, Bennett N, et al. Robust deep learning-based protein sequence design using ProteinMPNN[J]. Science, 2022, 378(6615): 49-56. That is: Dauparas J et al. Robust protein sequence design based on deep learning ProteinMPNN[J]. Science, 2022, 378(6615): 49-56), ESM-IF1 (HsuC, Verkuil R, LiuJ, et al. Learning inverse folding from millions of predicted structures[C]. International conference on machine learning. PMLR, 2022: 8946-8970. That is: Hsu C et al. Learning reverse folding from millions of predicted structures[C]. International conference on machine learning. PMLR, 2022: 8946-8970), CarbonDesign (RenM, Yu C, Bu D,et al.Accurate and robust protein sequence design with CarbonDesign[J].Nature Machine Intelligence,2024,6(5):536-547. That is: RenM et al. Accurate and robust protein sequence design with CarbonDesign[J].Nature Machine Intelligence,2024,6(5):536-547), ProGen (MadaniA,Krause B,Greene ER,etal.Large language models generate functional protein sequences across diverse families[J].Nature Biotechnology,2023,41(8):1099-1106. That is: MadaniA et al. Large language models generate functional protein sequences across diverse families[J].Nature Biotechnology,2023,41(8):1099-1106) ProteinMPNN uses a message passing mechanism graph neural network combined with an encoder-decoder architecture to restore protein sequences from the main chain skeleton.ESM-IF1 uses a geometric vector-aware graph neural network to learn the reverse folding pattern of proteins from the protein structures predicted by AlphaFold2. CarbonDesign uses multiple Inverseformer modules as the network backbone and designs different prediction tasks, including single and multiple amino acid loss, as well as side chain structure loss. ProGen uses protein language model technology to generate protein sequences of different families.

[0004] In summary, existing protein sequence design methods are still far from meeting the requirements of practical applications in terms of design accuracy and robustness, and are in urgent need of improvement. Summary of the Invention

[0005] The purpose of the present invention is to provide a protein sequence design method based on biomolecular interaction domain enhancement to overcome the defects of the prior art.

[0006] In order to achieve the above object, the technical solution adopted by the present invention is as follows:

[0007] The protein sequence design method based on biomolecular interaction domain enhancement includes the following steps:

[0008] S1. Input the three-dimensional coordinate information C of the protein backbone to be sequenced, which is of size L×N×3, where L is the number of residues, N is the number of backbone atoms of each residue, and 3 is the Cartesian space coordinate of each atom;

[0009] S2. Obtain protein sequences that contact biomolecules and interaction domains Where N is the total number of sequences in the sequence set D, seq i represents the i-th sequence in the sequence set D, domain i =[start,end] represents the domain interval of the i-th sequence in the sequence set D, where 1≤start≤end;

[0010] S3. Cluster the protein sequence D into M clusters and extract the representative sequences of all clusters to form a non-redundant sequence set and the corresponding domain interval set Where M is the total number of sequences in the non-redundant sequence set S;

[0011] S4, predicting that each protein sequence seq in the sequence set S is of length G i The three-dimensional structure information is represented as a PDB file, and the atomic coordinates of the three-dimensional structure information are marked as Coord;

[0012] S5. Calculate all the PDB files to obtain each protein sequence seq i The eight-state secondary structure SS={s1,s2,…,s j ,…,s G}, where s j ∈(H,B,E,G,I,T,S,L) represents protein sequence i The secondary structure state of the jth residue in ;

[0013] S6. Calculate all protein sequence seq i Each protein sequence in the PDB file i The solvent accessible surface area SA = {a1, a2, ..., a j ,…,a G}, where a is a non-negative real value j Protein seq i The solvent accessible surface area of ​​the jth residue of ;

[0014] S7. Obtain each protein sequence seq in the sequence set S. i Functional annotation information file in, represents the jth functional annotation information in FUN, and N represents the total number of functional annotation information in FUN;

[0015] S8. Use the low-rank adaptive algorithm LoRA to adjust the parameter weights of the last 10 layers of the transformer block of the universal multimodal protein language ESM3, and record the ESM3 model as the EiRA model;

[0016] S9, using the sequence set S as a training data set to train the EiRA model, using the sequence information seq i , three-dimensional structure information Coord, secondary structure information SS, solvent accessibility information SA and functional annotation information FUN as input features;

[0017] S10, using standard random interval mask prediction as a training task, when the first set ratio, the interaction domain interval domain i The residues in the interaction domain remain uncovered, and at the second set ratio, the interaction domain interval domain i The cross entropy loss of the masked residues in is multiplied by the weight coefficient, and the loss function is;

[0018]

[0019] Where, the random mask position n∈domain i, x is the token representation corresponding to the residue feature;

[0020] S11, using the AdamW algorithm to optimize the loss function, using pytorch's distributed data parallel algorithm to train the EiRA model, iterating on the training set S for a set number of cycles;

[0021] S12. Input the three-dimensional coordinate information C of the main chain skeleton of the protein to be sequenced into the trained EiRA model to obtain the output sequence of the target protein, which is the design result.

[0022] Furthermore, in step S2, all protein sequences and interaction domain regions that contact with biomolecules are downloaded from the UniProtKB database.

[0023] Furthermore, in step S3, the MMseqs2 tool with a clustering threshold of 0.4 is used to cluster the protein sequence D into M clusters and the representative sequences of all clusters are extracted.

[0024] Furthermore, in step S4, the AlphaFold2 tool is used to predict the protein sequence seq of length G in the sequence set S. i The three-dimensional structure information is represented as a PDB file.

[0025] Furthermore, in step S5, the DSSP tool is used to calculate all the PDB files to obtain each protein sequence seq i The eight-state secondary structure.

[0026] Furthermore, in step S6, the Shrake-Rupley algorithm is used to calculate the seq of each protein. i of solvent accessible surface area.

[0027] Furthermore, in step S7, the InterProScan program is used to search the Pfam, SFLD, PANTHER, PROSITE and CATH databases to obtain each protein sequence seq in the sequence set S. i Functional annotation information file.

[0028] Furthermore, the rank of the low-rank adaptive algorithm in step S8 is 8, the a parameter is 32, and the loss rate is 0.05.

[0029] Furthermore, the standard random interval mask prediction is used as a training task in step S10, specifically: a certain number of residues in a protein interval of length 7 are randomly masked and then predicted, and the total masked residue ratio in the sequence is 20%.

[0030] Further, the learning rate of the AdamW algorithm in the step S10 is 1e-4, the weight decay is 0.01, the gradient clipping is 1.0, beta1=0.9, and beta2=0.95.

[0031] Compared with the prior art, the advantages of the present application are that: on the one hand, the present application utilizes the multi-modal information of massive proteins, including sequence, structure and functional features, to provide a data basis for improving the accuracy of protein sequence design; on the other hand, more attention is explicitly given to protein functional domains, and the carefully designed learning task enables the model to generate more robust and reasonable functional protein sequences. BRIEF DESCRIPTION OF DRAWINGS

[0032] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or the prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor.

[0033] Figure 1 is the principle diagram of the protein sequence design method based on biomolecular interaction domain enhancement of the present application.

[0034] Figure 2 is the primary sequence diagram of the protein 1a73_A obtained by the present application. DETAILED DESCRIPTION

[0035] The preferred embodiments of the present application will be described in detail below with reference to the drawings, so that the advantages and features of the present application can be more easily understood by those skilled in the art, and the protection scope of the present application can be more clearly and definitely defined.

[0036] Referring to Figure 1 The present embodiment discloses a protein sequence design method based on biomolecular interaction domain enhancement, taking the sequence design of protein 1a73_A as an example, the present embodiment includes the following steps:

[0037] Step S1, input the protein backbone three-dimensional coordinate information C to be designed with a size of LxNx3, wherein L is the number of residues, N is the number of backbone atoms of each residue, and 3 is the Cartesian space coordinates of each atom.

[0038] Step S2, obtain the protein sequence and the interaction domain interval from the UniProtKB database which contacts with biomolecules (such as DNA, RNA and protein, etc.) irepresents the i-th sequence in the sequence set D, domain i =[start,end] represents the domain interval of the i-th sequence in the sequence set D, where 1≤start≤end.

[0039] Step S3: Use the MMseqs2 tool with a clustering threshold of 0.4 to cluster the protein sequence D into M clusters and extract the representative sequences of all clusters to form a non-redundant sequence set. and the corresponding domain interval set Where M is the total number of sequences in the non-redundant sequence set S.

[0040] Step S4: Use AlphaFold2 to predict the protein sequence seq of length G in the sequence set S. i The three-dimensional structure information is represented as a PDB file, and the atomic coordinates of the three-dimensional structure information are marked as Coord.

[0041] Step S5: Use the DSSP tool to calculate all PDB files to obtain each protein sequence seq i The eight-state secondary structure SS={s1,s2,…,s j ,…,s G}, where s j ∈(H,B,E,G,I,T,S,L) represents protein sequence i The secondary structure state of the jth residue in ;

[0042] Step S6: Calculate all protein sequences using the Shrake-Rupley algorithm i Each protein sequence in the PDB file i The solvent accessible surface area SA = {a1, a2, ..., a j ,…,a G}, where a is a non-negative real value j Protein seq i The solvent accessible surface area of ​​the jth residue in .

[0043] Step S7: Use the InterProScan program to search the Pfam, SFLD, PANTHER, PROSITE, and CATH databases to obtain each protein sequence seq in the sequence set S. i Functional annotation information file in, represents the jth functional annotation information in FUN, and N represents the total number of functional annotation information in FUN.

[0044] Step S8: Use the low-rank adaptive algorithm LoRA with a rank of 8, an a parameter of 32, and a dropout rate of 0.05 to adjust the parameter weights of the last 10 layers of the transformer block of the universal multimodal protein language ESM3, and record the ESM3 model as the EiRA model.

[0045] Step S9: Use the sequence set S as a training data set to train the EiRA model, and use the sequence information seq i , three-dimensional structure information Coord, secondary structure information SS, solvent accessibility information SA and functional annotation information FUN are used as input features.

[0046] Step S10, using standard random interval mask prediction as a training task (i.e., randomly masking a certain number of residues of a protein interval of length 7 and then making predictions, with the total masked residue ratio in the sequence being 20%), when the first set ratio is 60%, the interaction domain interval domain i The residues in the interaction domain remain uncovered. At the second setting ratio of 40%, the interaction domain interval domain i The cross entropy loss of the masked residues in is multiplied by the weight coefficient 1.4 to achieve the effect of enhanced learning of the domain. The loss function is;

[0047]

[0048] Where, the random mask position n∈domain i , x is the token representation corresponding to the residue feature.

[0049] Step S11, using the AdamW algorithm with a learning rate of 1e-4, a weight decay of 0.01, a gradient clipping of 1.0, β1=0.9, and β2=0.95 to optimize the loss function, using pytorch's distributed data parallel algorithm to train the EiRA model on 4 NVIDIA V10032G graphics cards, and iterating on the training set S for a set number of cycles (e.g., 15).

[0050] Step S12: Input the three-dimensional coordinate information C of the protein backbone to be sequenced into the trained EiRA model to obtain the sequence of the target protein. The sequence of the target protein 1a73_A is the design result. The primary sequence of the protein 1a73_A is as follows: Figure 2 shown.

[0051] The present invention provides a highly accurate and robust protein sequence design method based on enhanced biomolecular interaction domains. On the one hand, the present invention utilizes the multimodal information of massive proteins, including sequence, structure and functional characteristics, providing a data basis for improving the accuracy of protein sequence design; on the other hand, it explicitly pays more attention to the functional domains of proteins, and carefully designed learning tasks enable the model to generate more robust and reasonable functional protein sequences.

[0052] Although the embodiments of the present invention are described in conjunction with the accompanying drawings, the patent owner may make various changes or modifications within the scope of the appended claims. As long as they do not exceed the scope of protection described in the claims of the present invention, they should be within the scope of protection of the present invention.

Claims

1. A protein sequence design method based on biomolecular interaction domain enhancement, characterized in that: The following steps are involved: S1. Input the three-dimensional coordinate information C of the protein backbone to be sequenced, which is of size L×N×3, where L is the number of residues, N is the number of backbone atoms of each residue, and 3 is the Cartesian space coordinate of each atom; S2. Obtain protein sequences that contact biomolecules and interaction domains Where N is the total number of sequences in the sequence set D, seq i represents the i-th sequence in the sequence set D, domain i =[start,end] represents the domain interval of the i-th sequence in the sequence set D, where 1≤start≤end; S3. Cluster the protein sequence D into M clusters and extract the representative sequences of all clusters to form a non-redundant sequence set and the corresponding domain interval set Where M is the total number of sequences in the non-redundant sequence set S; S4, predicting that each protein sequence seq in the sequence set S is of length G i The three-dimensional structure information is represented as a PDB file, and the atomic coordinates of the three-dimensional structure information are marked as Coord; S5. Calculate all the PDB files to obtain each protein sequence seq i The eight-state secondary structure SS={s1,s2,…,s j ,…,s G }, where s j ∈(H,B,E,G,I,T,S,L) represents protein sequence i The secondary structure state of the jth residue in ; S6. Calculate all protein sequence seq i Each protein sequence in the PDB file i The solvent accessible surface area SA = {a1, a2, ..., a j ,…,a G }, where a is a non-negative real value j Protein seq i The solvent accessible surface area of ​​the jth residue of ; S7. Obtain each protein sequence seq in the sequence set S. i Functional annotation information file in, represents the jth functional annotation information in FUN, and N represents the total number of functional annotation information in FUN; S8. Use the low-rank adaptive algorithm LoRA to adjust the parameter weights of the last 10 layers of the transformer block of the universal multimodal protein language ESM3, and record the ESM3 model as the EiRA model; S9, using the sequence set S as a training data set to train the EiRA model, using the sequence information seq i , three-dimensional structure information Coord, secondary structure information SS, solvent accessibility information SA and functional annotation information FUN as input features; S10, using standard random interval mask prediction as a training task, when the first set ratio, the interaction domain interval domain i The residues in the interaction domain remain uncovered, and at the second set ratio, the interaction domain interval domain i The cross entropy loss of the masked residues in is multiplied by the weight coefficient, and the loss function is; Where, the random mask position n∈domain i , x is the token representation corresponding to the residue feature; S11, using the AdamW algorithm to optimize the loss function, using pytorch's distributed data parallel algorithm to train the EiRA model, iterating on the training set S for a set number of cycles; S12. Input the three-dimensional coordinate information C of the main chain skeleton of the protein to be sequenced into the trained EiRA model to obtain the output sequence of the target protein, which is the design result.

2. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In the step S2, all protein sequences and interaction domain regions that contact with biomolecules are downloaded from the UniProtKB database.

3. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S3, the protein sequence D is clustered into M clusters using the MMseqs2 tool with a clustering threshold of 0.4, and representative sequences of all clusters are extracted.

4. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S4, the AlphaFold2 tool is used to predict each protein sequence seq of length G in the sequence set S. i The three-dimensional structure information is represented as a PDB file.

5. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S5, the DSSP tool is used to calculate all the PDB files to obtain each protein sequence seq i The eight-state secondary structure.

6. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S6, the Shrake-Rupley algorithm is used to calculate the sequence of each protein i of solvent accessible surface area.

7. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S7, the InterProScan program is used to search the Pfam, SFLD, PANTHER, PROSITE and CATH databases to obtain each protein sequence seq in the sequence set S. i Functional annotation information file.

8. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: The rank of the low-rank adaptive algorithm in step S8 is 8, the a parameter is 32, and the loss rate is 0.

05.

9. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S10, standard random interval mask prediction is used as a training task, specifically: a certain number of residues in a protein interval of length 7 are randomly masked and then prediction is performed, and the total masked residue ratio in the sequence is 20%.

10. The protein sequence design method based on biomolecular interaction domain enhancement according to claim 1, characterized in that: In step S10, the learning rate of the AdamW algorithm is 1e-4, the weight decay is 0.01, the gradient clipping is 1.0, β1=0.9, and β2=0.95.

Images