Protein sequence design method based on biomolecule interaction structural domain enhancement

Abstract

The invention discloses a protein sequence design method based on biomolecular interaction structural domain enhancement. The method comprises the following steps: firstly, inputting three-dimensional coordinate information of a protein main chain skeleton to be subjected to sequence design, wherein the size of the protein main chain skeleton is L * N * 3; obtaining a protein sequence in contact with the biomolecules and an interaction structural domain interval; clustering the obtained sequences and taking out a representative sequence of each cluster as a training set; extracting a three-dimensional structure, a secondary structure, solvent accessibility and function annotation feature representation of each training sample; fine-tuning the last ten layers of transformer modules of the universal multi-modal protein language model ESM3 by using an LoRA algorithm, and giving a larger weight for the loss of mask residues located in an interaction structural domain interval; and inputting atomic coordinates of a protein main chain skeleton to be subjected to sequence design into the trained model to obtain a target sequence. On one hand, multi-modal information of massive proteins is utilized; on the other hand, a more robust and reasonable functional protein sequence can be generated.

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 biomolecular interaction domain enhancement. Background Art

[0002] Protein sequence design refers to designing a primary amino acid sequence that can fit the topological structure based on the protein backbone. Since artificially designed proteins have artificially specified functions, they have a broad and far-reaching impact on the field of life sciences. For example, by designing specific protein sequences, new therapeutic drugs can be developed or the efficacy of existing drugs can be enhanced; by designing more stable and efficient enzymes, the efficiency of industrial processes can be improved and costs can be reduced; it can be used to design antigens with strong immunogenicity that can induce effective immune responses. Therefore, the precise design of protein sequences is of great significance to life science fields such as drug development, vaccine development, and enzyme engineering. With the rise of artificial intelligence technology, the use of deep learning to design new protein sequences from scratch has received 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: Ren M et al. Accurate and robust protein sequence design with CarbonDesign[J]. Nature Machine Intelligence, 2024, 6(5): 536-547), ProGen (Madani A, Krause B, Greene ER, et al. Large language models generate functional protein sequences across diverse families[J]. Nature Biotechnology, 2023, 41(8): 1099-1106. That is: Madani A et al. Large language models generate functional protein sequences across diverse families[J]. Nature Biotechnology, 2023, 41(8): 1099-1106) Protein MPNN 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 mode of proteins from the protein structure predicted by AlphaFold2. CarbonDesign uses multiple Inverseformer modules as the network backbone and designs different prediction tasks, including single and multiple amino acid losses, and side chain structure losses. ProGen uses protein language model technology to generate protein sequences of different families.

[0004] In summary, the existing protein sequence design methods are still far from 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 main chain skeleton of the protein to be sequenced, which is L×N×3 in size, where L is the number of residues, N is the number of main chain atoms of each residue, and 3 is the Cartesian space coordinate of each atom;

[0009] S2. Obtain protein sequences that contact biological molecules 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, clustering the protein sequence D into M clusters and extracting 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 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;

[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 seq i The secondary structure state of the jth residue in ;

[0013] S6. Calculate all protein sequences 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, obtaining 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, using 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 the ESM3 model is recorded 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, at a 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, and iterating on the training set S for a set of cycles;

[0021] S12, inputting 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 biological molecules are specifically 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 extract representative sequences of all clusters.

[0024] Furthermore, in step S4, the AlphaFold2 tool is used to predict that each 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 sequence 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 step S10 uses standard random interval mask prediction as a training task, specifically: randomly masking a certain number of residues in a protein interval of length 7 and then performing prediction, and the total masked residue ratio in the sequence is 20%.

[0030] Furthermore, 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.

[0031] Compared with the prior art, the advantages of the present invention are: 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 the carefully designed learning tasks enable the model to generate more robust and reasonable functional protein sequences. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0033] Figure 1 It is a schematic diagram of the protein sequence design method based on biomolecular interaction domain enhancement of the present invention.

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

[0035] The preferred embodiments of the present invention are described in detail below in conjunction with the accompanying drawings so that the advantages and features of the present invention can be more easily understood by those skilled in the art, thereby making a clearer and more definite definition of the protection scope of the present invention.

[0036] See also Figure 1 As shown, this embodiment discloses a protein sequence design method based on biomolecular interaction domain enhancement. Taking the sequence design of protein 1a73_A as an example, this embodiment includes the following steps:

[0037] Step S1, inputting the three-dimensional coordinate information C of the main chain skeleton of the protein to be sequenced with a size of L×N×3, wherein L is the number of residues, N is the number of main chain atoms of each residue, and 3 is the Cartesian space coordinate of each atom.

[0038] Step S2: Obtain protein sequences that contact biological molecules (such as DNA, RNA, and proteins) from the UniProtKB database and interaction domains Where N is the total number of sequences in the sequence set D, seq irepresents the i-th sequence in the sequence set D, domain i =[start,end] represents the structural 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 seq i The secondary structure state of the jth residue in ;

[0042] Step S6: Calculate all protein sequences seq 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, a parameter of 32, and a loss rate of 0.05 to adjust the parameter weights of the last 10 layers of transformer blocks 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 in 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, β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 (e.g., 15) cycles.

[0050] Step 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 sequence of the target protein as output. The sequence of the target protein 1a73_A is the design result. The primary sequence of the obtained 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 biomolecular interaction domain enhancement. On the one hand, the present invention utilizes the multimodal information of massive proteins, including sequence, structure and functional characteristics, to provide a data basis for improving the accuracy of protein sequence design; on the other hand, it explicitly pays more attention to protein functional domains, 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 variations or modifications within the scope of the appended claims. As long as they do not exceed the protection scope described in the claims of the present invention, they should be within the protection scope 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 main chain skeleton of the protein to be sequenced, which is L×N×3 in size, where L is the number of residues, N is the number of main chain atoms of each residue, and 3 is the Cartesian space coordinate of each atom; S2. Obtain protein sequences that contact biological molecules 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, clustering the protein sequence D into M clusters and extracting 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 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; 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 seq i The secondary structure state of the jth residue in ; S6. Calculate all protein sequences 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, obtaining 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, using 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 the ESM3 model is recorded 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, at a 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, and iterating on the training set S for a set of cycles; S12, inputting 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 biological molecules 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 MMseqs2 tool with a clustering threshold of 0.4 is used to cluster the protein sequence D into M clusters and extract representative sequences of all clusters.

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 that 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 predicted, 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