A method for predicting heterodimeric interchain residue contacts

Abstract

The application discloses a kind of heterodimer interchain residue contact prediction method.The present application is aimed at the deficiency of existing method in prediction accuracy, long-range dependence modeling capability and generalization, and proposes a kind of deep neural network integrating multiple features and attention mechanism.Specifically, the application integrates multiple features such as protein language model as network model input, then the network adopts efficient channel attention (ECA) and spatial attention (SA) module and KAN convolution network module, effectively captures the local and global dependence features of heterodimer, so as to predict the interchain residue contact of heterodimer.Experiments show that the prediction accuracy of the application on the benchmark dataset is significantly better than that of existing methods, and the model has high robustness.The application can be widely used in the field of protein heterodimer interchain residue contact prediction and protein structure prediction, further promoting protein function research and protein drug development.

Description

TECHNICAL FIELD

[0001] The present application relates to protein-protein interaction, and in particular to a method for predicting inter-chain residue contacts of heterodimers. BACKGROUND

[0002] Protein heterodimers are complexes composed of two different types of protein monomers (different lengths, different structures), and heterodimers are usually composed of protein chains with different functional and structural properties. These different protein monomers are often complementary in function, and the interaction between the two is crucial for achieving specific biological functions. For example, many key biological molecules such as receptors, enzymes and transcription factors exist in the form of heterodimers, which regulate cell signal transduction, metabolic processes and gene expression by forming heterodimer complexes. The formation of heterodimers is usually driven by protein-protein interactions (PPIs), and these interactions often have high specificity and selectivity. The interaction between protein chains not only depends on electrostatic forces, hydrogen bonds and hydrophobic interactions, but also involves the fine docking of molecular structures. Therefore, understanding and predicting the inter-chain residue contacts of heterodimers is of great significance for understanding the structure and function of protein complexes. Compared with homodimers, the inter-chain contact prediction of heterodimers is more difficult and faces greater challenges, because there is more heterogeneity between the two chains of heterodimers, and the differences in structure and functional properties make the interaction pattern more complex.

[0003] From the currently published methods for predicting inter-chain residue contacts of heterodimers, the overall prediction accuracy is low, and the accuracy is less than 50% on most prediction samples. The reasons are as follows: first, it is difficult to build high-quality joint MSAs, especially when the sequence difference between heterodimers in different species is large, it is often difficult to find enough homologous sequences, and when the quality of MSAs is poor and the sequence pair matching is inaccurate, the prediction results will be affected by noise. Second, the existing models are not well integrated, different models are good at capturing different features, but lack of integration of existing models. In general, although the existing methods for predicting inter-chain residue contacts of heterodimers have made some progress, they still face problems such as MSA pairing, noise processing, model integration, and low overall prediction accuracy. SUMMARY

[0004] This invention provides an AttCON-Hetero method for predicting inter-chain residue contacts in heterodimers, based on deep learning and a protein language model. The method employs multiple feature fusion techniques to process input features. The network model combines channel attention (ECA) and spatial attention (SA) modules with a KAN convolutional network module, achieving higher performance in predicting inter-chain residue contacts and providing highly valuable information for dimer structure prediction, functional prediction, and drug design.

[0005] The technical solution of this invention is: a method for predicting interchain residue contacts in heterodimers, comprising the following steps:

[0006] S1. Screening out heterodimer structure information (Heterodimer-PDB) from protein databases;

[0007] S2. Divide the data obtained in step S1 into training set, validation set, and test set. Based on the heterodimer structure information, preprocess the Heterodimer-PDB structure data to extract the structural information MonomerA-PDB / MonomerB-PDB and sequence information MonomerA–FATSA / MonomerB-FASTA of the two monomers of the heterodimer, with sequence lengths L respectively. A and L B And directly spliced ​​together to form a line of length L A +L B The combined sequence FASTA; based on the heterodimer structure information, Heterodimer-PDB was used to calculate the Euclidean distance between the two monomeric amino acids of the heterodimer, thus obtaining the inter-chain contact between the heterodimers (matrix size L). A ×L B );

[0008] S3. Using the two monomer sequence information MonomerA–FATSA / MonomerB-FASTA generated in step S2, the HHblits multiple sequence alignment search algorithm is used to search the Uniclust30 sequence database to obtain the multiple sequence alignment information (MSAs) of these two monomers respectively. The MSAs of the two monomers are matched using the TaxID species classification number to ensure that the selected homologous sequences come from the same species, and finally a joint multiple sequence alignment Paired MSAs is generated. This method can integrate information from two monomers and capture the co-evolutionary signals between them.

[0009] S4, using the paired MSA information generated in step S3, using ESM-MSA and ESM2 two protein language models to extract line attention features respectively; using CCMpred algorithm to extract coevolution information features in paired MSAs; using the structure information of the two monomers of the heterodimer in step S2 to calculate the distance map of the monomer amino acid chain respectively; using the joint sequence FASTA obtained in step S2, using the PSI-BLAST tool to generate the position-specific scoring matrix (PSSM);

[0010] S5, the 7 types of features generated in step S4 are fused to generate a fusion feature matrix, which improves the diversity of input features;

[0011] S6, input the features processed in step S5 into the model AttCON-Hetero composed of deep residual network ResNet and KAN (Kolmogorov-Arnold Networks) convolutional network, said model uses efficient channel attention ECA and spatial attention SA module to capture local and global dependent features of the feature matrix, and uses KAN convolution with self-adaptive convolution kernel weight to output the predicted contact probability matrix, and finally outputs the probability matrix of amino acid contact between the two chains of the heterodimer protein (matrix size L A ×L B );

[0012] S7, based on the network output of step S6, use the contact matrix Inter-ChainContact between the heterodimers in step S2 as the label for training, then use the validation set of S2 to verify the performance of the model obtained by each batch training, and save the model parameters with the best indicators obtained on the validation set as the weight of the final classification model. Finally, verify the performance of the saved model on the test set of S2.

[0013] S1 selects the heterodimer from the PDB database according to the following conditions: (1) the aggregation state should be heterodimer; (2) the biological assembly contains only two single chains; (3) the resolution of X-ray structure should be better than (4) the length of the single chain should be greater than 50, and the total length of the complex should be less than 500. (5) using 40% sequence similarity threshold to cluster the heterodimer, and selecting the best resolution heterodimer from each cluster.

[0014] As a further description of the above scheme: S2 extracts the structure information MonomerA-PDB / MonomerB-PDB and sequence information of the heterodimer monomer, and calculates the contact matrix between the amino acids of the two monomers of the heterodimer, including the following steps:

[0015] S2.1, extract the information of heterodimer monomer and sequence information:

[0016] Extract the structure information and sequence information of the two chains of the heterodimer: for the qualified dimer, extract the three-dimensional structure information of the two chains from its PDB file, save them as Monomer-PDB A and Monomer-PDB B respectively, extract the amino acid sequence information MonomerA-FATSA and MonomerB-FASTA, and combine MonomerA-FATSA and MonomerB-FASTA into a sequence FASTA; during the extraction process, remove the residues with amino acid label UNK;

[0017] S2.2, calculate the Euclidean distance and contact matrix between the amino acid residues of the two chains of the heterodimer:

[0018] For residues i and j (i from monomer A; j from monomer B), two types of distance calculation rules are defined:

[0019]

[0020] i∈MonomerA,j∈MonomerB (1)

[0021] In formula (1), d is the nearest heavy atom distance, where r a and r b represent the three-dimensional coordinates of the heavy atoms in residues i and j, respectively, and ||r a -r b is the Euclidean distance of the nearest heavy atoms between residues i and residues j.

[0022]

[0023] In formula (2), d is the distance between the beta carbon atoms (C β -distance) of residues i and residues j, is the C β atomic coordinate of residue i; according to formula (1) and formula (2), two types of contact matrices are calculated, if or is less than or equal to , the value of the corresponding position of the contact matrix is 1, otherwise 0, and they are saved as Inter-Chain Contact(Heavy Atom) and Inter-Chain Contact(C β ), respectively. The matrix size is L A ×L BThe above matrices are used as labels of the model for comparison of the experimental results.

[0024] As a further description of the above scheme: S3 specifically includes the following steps:

[0025] S3.1, using HHblits to search the Uniclust30 database to obtain MSAs:

[0026] Using the monomer sequence information generated in step S2.1, the Uniclust30 database is preliminarily searched using the HH-suite 3.30 tool;

[0027] S3.2, matching the MSAs of the two monomers to generate a joint multiple sequence alignment (Paired MSAs):

[0028] The MSAs of monomer A and monomer B are matched using the species classification number TaxID of each monomer; it is ensured that the two monomers come from the same biological species or related homologous species; according to the TaxID information, the matching results of monomer A and monomer B are integrated to generate a joint multiple sequence alignment Paired MSAs; this joint MSA contains homologous sequence information from the two monomers, providing high-quality input data for subsequent feature extraction.

[0029] As a further description of the above scheme: S4 uses the joint multiple sequence alignment information (Paired MSAs) and other feature sources to extract various features, and the detailed steps are as follows:

[0030] S4.1, extract the row attention feature (F ESM-MSA and F ESM2 )

[0031] F ESM-MSA , i.e. ESM-MSA row attention feature: based on the joint multiple sequence alignment information (Paired MSAs) generated in step S3, the ESM-MSA model (the model used in this application is esm_msa1_t12_100M_UR50S) is used to extract the row attention feature in the sequence. The ESM-MSA model is specifically optimized for multiple sequence alignment, which can capture the co-evolution signals in MSAs and the differences between sequences, and generate a high-quality row attention feature matrix (ESM-MSA Row Attention).

[0032] F ESM2ESM2 row attention feature: Extract row attention feature by ESM2 model (the model used in the present application is esm2_t6_8M_UR50D) using joint sequence FASTA. ESM2 model performs well in capturing long-range dependencies and local properties of protein sequences, and the generated row attention feature matrix ESM2 RowAttention has strong generalization and prediction accuracy.

[0033] S4.2, extract co-evolution information feature (F CCMpred , F DCA-APC and F DCA-DI )

[0034] F CCMpred CCMpred co-evolution feature: Based on Paired MSAs, use CCMpred algorithm to calculate the probability distribution of inter-sequence co-evolution, and generate a diagonal symmetric co-evolution matrix;

[0035] F DCA-APC and F DCA-DI : Further based on the Direct Coupling Analysis (DCA) method, respectively generate DCA feature matrix with average product correction (DCA-APC) and uncorrected DCA feature matrix (DCA-DI), these features reveal the strength and pattern of sequence co-evolution, which are used to describe the interaction between residues;

[0036]

[0037] F DCA-DI = DCA (4)

[0038] DCA in formula (3) is the result of direct coupling analysis, Row Mean, Col Mean and Overall Mean are the row mean, column mean and average of all elements of the DCA matrix respectively;

[0039] S4.3, extract two monomer intra-chain distance features (F intra )

[0040] Based on the two monomer structure information (MonomerA-PDB / MonomerB-PDB) generated in step S2, the Euclidean distance matrix (Intra-Chain Distance) between the amino acids in the two monomer chains is calculated respectively; for each pair of residues (a, b) in the monomer chain, the distance can be expressed as:

[0041]

[0042] In formula (5), (x a , ya a ) and (x b ,y b ,z b ) are the three-dimensional coordinates of the a-th and b-th residues, respectively; the intra-chain distance matrix generated by monomer A is D A , and the intra-chain distance matrix generated by monomer B is D B , which can be expressed as follows:

[0043]

[0044] S4.4, generating a position-specific scoring matrix (F PSSM )

[0045] Using the joint sequence information (FASTA) in step S2, the position-specific scoring matrix (PSSM) of the sequence is calculated by the PSI-BLAST tool; the PSSM captures the substitution patterns and conservation information of amino acids by multiple iterative searches of the UniProt database, and the generated matrix can be used for feature enhancement and model optimization.

[0046]

[0047] In formula (7), f ij is the frequency of amino acid j at position i, b j is the background frequency.

[0048] As a further description of the above scheme: the detailed steps of S5 for generating the fusion feature are as follows:

[0049] The seven features generated in step S4 are concatenated according to the channel dimension to generate the fusion feature F contact ;

[0050] F contact = [F ESM-MSA , F ESM2 , F CCMpred , F DCA-APC , F DCA-DI , F intra , F PSSM ] (8)

[0051] In formula (8), the dimension of F contact is (L A + L B ) × (L A + L B ) × C, where L A and L B are the sequence lengths of monomer A and monomer B, respectively, and C represents the number of channels after fusion.

[0052] ​As a further description of the above scheme: S6 network model construction based on deep learning and protein language model includes the following steps:

[0053] The enhanced feature matrix F contact generated in step S5 is input into an AttCON-Hetero model based on deep residual network (ResNet) and KAN (Kolmogorov-Arnold Networks) convolution for feature learning, and the dimension of the feature matrix is (L A +L B )×(L A +L B )×C; the input feature matrix is standardized using instance normalization (Instance Normalization) to reduce the distribution deviation of the input features; C feature channels are compressed by two-dimensional convolution (2D Convolution); and the features of the channels are further extracted by the Maxout layer to obtain X input ; then X input is input into a residual module, each residual module is composed of the following steps, which effectively extracts local and global features by step-by-step standardization, convolution, and attention module, while maintaining the consistency of the features;

[0054] First, the input feature X input is standardized to unify the feature scale:

[0055]

[0056] where Norm is the instance normalization operation, and then two-dimensional convolution is used to extract preliminary features:

[0057]

[0058] where Conv2D is a two-dimensional convolution operation, in order to capture multi-scale feature information, three different size convolution kernels are then used to further extract features, and the three groups of features are spliced to fuse different scale information:

[0059]

[0060] X multi =Contact(X1,X2,X3) (14)

[0061] where X1, X2, X3 are the results of the features extracted by different convolution kernels, Contact is a concatenation operation, and then X multi is standardized and two-dimensional convolution compressed according to formulas (9), (10) to obtain Subsequently, Efficient Channel Attention (ECA) and Spatial Attention (SA) are used for feature extraction:

[0062]

[0063] z=σ(Conv1D(y)) (16)

[0064]

[0065] f sa =Conv2D([AvgPool(X ECA ).;.MaxPool(X ECA (18)

[0066] X SA =X ECA ⊙σ(f sa (19)

[0067] Formulas (15) to (17) represent the processing flow of the channel attention mechanism: GAP is global max pooling, which averages the spatial dimensions of each channel to generate a global feature representation; Conv1D is a one-dimensional convolution operation; σ is the Sigmoid activation function; ⊙ represents element-wise multiplication. Formulas (18) to (19) represent the spatial attention mechanism: AvgPool and MaxPool are average pooling and max pooling operations, respectively; [·; ·] represents the concatenation of pooling results; X SA This is the result after spatial attention (SA) extracts features.

[0068]

[0069] Formula (20) adds the output of the attention mechanism to the input features through a residual connection, thereby enhancing gradient flow and avoiding the gradient vanishing problem;

[0070] Finally, there's the model's output layer, which takes X... output Standardize according to formula (9) to obtain X output_norm The output is then produced through two layers of KAN convolution. KAN convolution ConvKAN adaptively adjusts the weights of the convolution kernel, enabling the network to dynamically adjust its learning process according to the features of the input data. This has significant advantages for capturing the heterogeneity and complex contact patterns between protein chains.

[0071] X KAN =ConvKAN(X output_norm ) (twenty one)

[0072] ConvKAN is KAN convolution, X KAN is the result after KAN convolution processing. The output of formula (21) uses a KAN convolution layer and a Sigmoid classifier, and outputs a probability matrix (contact or not), and the output matrix dimension is L A x L B , which represents the contact probability P contact of each pair of residues:

[0073] P contact = Sigmoid(ConvKAN(X KAN )) (22)

[0074] As a further description of the above scheme: the S7 model training process, model performance evaluation and model test steps are as follows:

[0075] S7.1 Use of loss function in model training and optimization steps of optimizer

[0076] The contact prediction matrix output in step S6 is taken as the output of the network, and the inter-chain residue contact matrix (Inter-Chain Contact) obtained in step S2 is taken as the label label (in the experiment, Inter-Chain Contact (C β ) is preferred as the label, and if the Inter-Chain Contact (C β ) part sample calculation is problematic, Inter-Chain Contact (HeavyAtom) is uniformly selected as the label), and the model is trained; the training process is optimized by Focal Loss loss function, using AdamW optimizer combined with cosine annealing learning rate scheduling;

[0077] S7.2 Performance evaluation and model selection

[0078] After each batch of training, the model is evaluated for performance using the validation set provided in step S2, and the model performance is evaluated by calculating the Top n accuracy, Top-L / k accuracy and weighted evaluation index W score , and finally the model parameters with the best performance on the validation set are selected.

[0079] The application has the following characteristics: because there is more heterogeneity between the two chains of the heterodimer, the interaction mechanism is complex, and the existing method still has a large space for improvement in feature integration, model optimization, etc., and overall, the current method has low accuracy in predicting inter-chain residue contact prediction of heterodimer. The application effectively integrates and fuses traditional co-evolution features, monomer structure information and transformer row attention features generated by protein language models (ESM-MSA and ESM2), and optimizes the prediction model, uses Efficient Channel Attention (ECA) and Spatial Attention (SA) attention modules to efficiently extract features, wherein ECA is used to weight channel features, SA is used to capture spatial information and integrate channel features, and finally uses the KAN convolutional network module to adaptively adjust the convolution kernel weight capability, efficiently output the inter-chain residue contact probability map. Compared with the existing prediction method, more accurate inter-chain residue contact information can be obtained.

[0080] Compared with the prior art, the application has the following technical effects: the application integrates the attention features generated by the protein language model (ESM-MSA and ESM2), which enhances the diversity of input data and the robustness of the model compared with traditional protein structure features (such as distance map, co-evolution feature, etc.). Subsequently, the network adopts the efficient channel attention (ECA) and spatial attention (SA) modules and the KAN convolutional network module to effectively capture the local and global dependence features of the heterodimer. Experiments show that the prediction accuracy of the application on the benchmark dataset is significantly better than that of the existing methods (such as DeepComplex, GLINTER, CDPred, DeepInter, etc.), and higher prediction performance can be achieved. The application can be widely used in the fields of protein heterodimer inter-chain residue contact prediction and protein structure prediction, and further promotes the research of protein function and the development of protein drugs. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1 is the flow chart of the method AttCON-Hetero inter-chain residue contact prediction of heterodimer;

[0082] Figure 2 is the residual module data processing flow chart of the method AttCON-Hetero;

[0083] Figure 3 is the attention module data processing flow chart of the method AttCON-Hetero;

[0084] Figure 4 is the PDB four-level structure diagram of two heterodimer samples in the training set;

[0085] Figure 5is the inter-chain residue contact prediction map of the test set sample without using the prediction method. DETAILED DESCRIPTION

[0086] The application will be further described in conjunction with the accompanying drawings and examples, but the content of the application is not limited to the scope described.

[0087] Example 1

[0088] An inter-chain residue contact prediction method of a heterodimer, the flow chart is as shown in Figure 1 The method specifically comprises the following steps:

[0089] S1, obtaining the heterodimer data from the PDB database, the schematic diagram of the three-dimensional structure of the heterodimer is as shown in Figure 4 The a and b in the figure are two different heterodimer samples, for example, the heterodimer in b is composed of two monomers in green and brown, and the inter-chain residue contact to be predicted is the contact information between the two monomers;

[0090] S2, dividing the heterodimer data meeting the conditions from the PDB database into a training set, a validation set and a test set, Table 1 shows the number of samples in the training set, the validation set and the test set, according to the information in these sets, obtaining the two monomer PDB structure (MonomerA-PDB / MonomerB-PDB) and sequence information (MonomerA-FASTA / MonomerB-FASTA / FASTA) of each sample and the inter-chain residue contact matrix (Inter-ChainContact) of each sample;

[0091] Table 1: Division of the data set used in the experiment

[0092]

[0093]

[0094] S3, using the monomer sequence information of all samples generated in step S2, searching the Uniclust30 sequence database using the HHblits multiple sequence alignment search algorithm to obtain the multiple sequence alignment information (MSAs) of all monomers, and using TaxID (Species Classification Number) to match the MSAs of each heterodimer sample, ensuring that the selected homologous sequences come from the same species, and pairing according to the species information, finally generating a joint multiple sequence alignment (Paired MSAs) for each sample.

[0095] S4, using the paired multiple sequence alignment information (Paired MSAs) of each heterodimer sample generated in step S3, using ESM-MSA and ESM2 two protein language models (PLMs) to extract line attention features (ESM-MSAAttention, ESM2-Attention) respectively; using the CCMpred algorithm to extract the coevolution information features (CCMpred, DCA-APC, DCA-DI) in the paired multiple sequence alignment information (Paired MSAs); using the structure information (MonomerA-PDB / MonomerB-PDB) of the two monomers of the heterodimer in step S2 to calculate the distance map (Intra-Chain DistanceA and Intra-Chain DistanceB) of the amino acid chain of the monomer respectively; using the joint sequence information (FATSA) obtained in step S2, using the PSI-BLAST tool to generate the position-specific scoring matrix (PSSM);

[0096] S5, the 7 types of features of each sample generated in step S4 are fused to generate a fusion feature matrix, and the diversity of the input features is improved;

[0097] S6, the features of the training samples processed in step S5 are input into a model (AttCON-Hetero) composed of a deep residual network (ResNet, Figure 2 and a KAN convolutional network, the model uses efficient channel attention (ECA) and spatial attention (SA) modules Figure 3 to capture local and global dependent features of the feature matrix, and uses a KAN convolution with self-adaptive adjustment of the convolution kernel weight to output a predicted contact probability matrix, and finally outputs a result as a probability matrix of the contact between two amino acids between the two chains of the heterodimer protein (the matrix size is L A ×L B ).

[0098] S7, based on the network output of step S6, the inter-chain residue contact prediction probability matrix is obtained, the contact matrix (Inter-Chain Contact) between the heterodimers of the training samples obtained in step S2 is used as a label for training, then the model performance of each batch of training is verified using the validation set of S2, and the model parameters with the best indicators obtained on the validation set are saved as the weight of the final classification model. Finally, the performance of the saved model is verified on the test set of S2.

[0099] Further, for S2, extract the training set, validation set and test set heterodimer monomer structure information (MonomerA-PDB / MonomerB-PDB) and sequence information, calculate the Euclidean distance between the amino acids of the two monomers of the heterodimer, and obtain the contact label data (Inter-Chain Contact).

[0100] S2.1, extract the three-dimensional structure information of the two chains from the PDB files of the training set, validation set and test set, save them as Monomer-PDB A and Monomer-PDB B respectively, extract the amino acid sequence information MonomerA-FATSA and MonomerB-FASTA, and combine MonomerA-FATSA and MonomerB-FASTA into a sequence FASTA. During the extraction process, remove the residues with amino acid label UNK.

[0101] S2.2, use the heterodimer PDB files of the training set, validation set and test set to calculate the Euclidean distance and contact matrix between the amino acid residues of the heterodimer chains:

[0102] For residues i and j (i from monomer A; j from monomer B), two types of distance calculation rules are defined:

[0103]

[0104] i∈MonomerA,j∈MonomerB(1)

[0105] In formula (1), d is the nearest heavy atom distance, where r a and r b represent the three-dimensional coordinates of the heavy atoms in residues i and j, respectively, and ||r a -r b is the Euclidean distance of the nearest heavy atoms between residues i and residues j.

[0106]

[0107] In formula (2), d is the distance between the beta carbon atoms (C β -distance) of residues i and residues j, is the C β atom coordinate of residue i; two types of contact matrices are calculated according to formula (1) and formula (2) respectively, and if or is less than or equal to The value of the contact matrix at the corresponding position is 1, otherwise 0, which are saved as Inter-Chain Contact (Heavy Atom) and Inter-Chain Contact (Light Atom) respectively. β The matrix size is L A x L B The above matrices are used as labels of the model for comparison of experimental results.

[0108] Further, the HHblits multiple sequence alignment search algorithm is used to search the Uniclust30 sequence database to obtain multiple sequence alignment information (MSAs) of two monomers of each sample in each dataset for S3, and TaxID (Species Classification Number) is used to match the MSAs of the two monomers to generate paired MSAs, including the following steps:

[0109] S3.1, using HHblits to search the Uniclust30 database to obtain MSAs:

[0110] Using the monomer sequence information generated in step S2.1, the Uniclust30 database is preliminarily searched using the HH-suite 3.30 tool. Uniclust30 is a database composed of UniProtKB sequences, all sequences are clustered based on 30% pairwise similarity, which can capture rich co-evolution signals. By setting appropriate search parameters (in this application, -diffinf-id 99-cov 50-n 3 is used), MSAs information of the two monomers is generated respectively.

[0111] S3.2, matching the MSAs of the two monomers to generate paired MSAs:

[0112] The MSAs of monomer A and monomer B are matched using the species classification number (TaxID) of each monomer. TaxID is used to identify the species to which each monomer belongs, and through the species information, it is ensured that the sequences of the two come from the same biological species or related homologous species. According to the TaxID information, the matching results of monomer A and monomer B are integrated to generate a paired MSA. This joint MSA contains homologous sequence information from the two monomers, providing high-quality input data for subsequent feature extraction.

[0113] Further, the detailed steps of extracting multiple features using paired multiple sequence alignment information (Paired MSAs) and other feature sources for S4 are as follows:

[0114] S4.1, extract row attention features (F ESM-MSA and FESM2 )

[0115] F ESM-MSA ESM-MSA row attention feature: Based on the paired MSAs generated in step S3, the ESM-MSA-1 model (protein language model, esm_msa1_t12_100M_UR50S is used in this application) is used to extract the row attention features in the sequence. The ESM-MSA model is specifically optimized for multiple sequence alignment, and by capturing the co-evolution signals and differences between sequences in MSAs, a high-quality row attention feature matrix (ESM-MSA Row Attention) is generated.

[0116] F ESM2 ESM2 row attention feature: Based on the joint sequence (FASTA), the row attention feature is extracted by the ESM2 model (protein language model, esm2_t6_8M_UR50D is used in this application). The ESM2 model performs well in capturing the long-range dependencies and local characteristics of protein sequences, and the generated row attention feature matrix (ESM2 Row Attention) has strong generalization and prediction accuracy.

[0117] S4.2, extract co-evolution information features (F CCMpred , F DCA-APC and F DCA-DI )

[0118] F CCMpred CCMpred co-evolution feature: Based on Paired MSAs, the CCMpred algorithm is used to calculate the probability distribution of co-evolution between sequences, and a diagonal symmetric co-evolution matrix is generated.

[0119] F DCA-APC and F DCA-DI : Further based on the Direct Coupling Analysis (DCA) method, the CCMpred algorithm is used to generate the DCA feature matrix (DCA-APC) and the uncorrected DCA feature matrix (DCA-DI) respectively. These features reveal the strength and pattern of sequence co-evolution, which are used to describe the interaction between residues.

[0120]

[0121] F DCA-DI = DCA (4)

[0122] DCA in formula (3) is the direct coupling analysis result, Row Mean, Col Mean and Overall Mean are the row mean, column mean and the average of all elements of the DCA matrix respectively.

[0123] S4.3, extract two monomer intra-chain distance features (F intra )

[0124] Based on the two monomer structure information (MonomerA-PDB / MonomerB-PDB) generated in step S2, the Euclidean distance matrix between the amino acids in the intra-chain of the two monomers is calculated respectively (Intra-Chain Distance). For each pair of residues (a, b) in the intra-chain of the monomer, the distance can be represented as:

[0125]

[0126] In formula (5), (x a ,y a ,z a ) and (x b ,y b ,z b ) are the three-dimensional coordinates of the a-th and b-th residues respectively. The intra-chain distance matrix generated by monomer A is D A , and the intra-chain distance matrix generated by monomer B is D B , which can be represented as follows:

[0127]

[0128] S4.4, generate position-specific scoring matrix (F PSSM

[0129] Using the joint sequence information (FASTA) in step S2, the position-specific scoring matrix (PSSM) of the sequence is calculated by the PSI-BLAST tool. PSSM captures the substitution patterns and conservation information of amino acids by multiple iterations of searching the UniProt database, and the generated matrix can be used for feature enhancement and model optimization.

[0130]

[0131] In formula (7), f ij is the frequency of amino acid j at position i, and b j is the background frequency.

[0132] Further, the detailed steps for generating the fusion feature S5 are as follows:

[0133] The seven features generated in step S4 are concatenated according to the channel dimension to generate the fusion feature F contact .​

[0134] F contact = [F ESM-MSA , F ESM2 , F CCMpred , F DCA-APC , F DCA-DI , F intra , F PSSM ] (8)

[0135] The dimension of F contact in formula (8) is (L A +L B )×(L A +L B )×C, wherein L A and L B are the sequence lengths of monomer A and monomer B, respectively, and C represents the number of channels after fusion, which is 308 in the present application.

[0136] Further, the network model based on deep learning and protein language model of S6 includes the following steps:

[0137] The feature matrix F contact generated in step S5 is input into an AttCON-Hetero model based on a deep residual network (ResNet) and KAN convolution for feature learning, and the dimension of the feature matrix is (L A +L B )×(L A +L B )×C. The input feature matrix is standardized using instance normalization (Instance Normalization) to reduce the distribution deviation of the input features. The feature channels are compressed from 308 channels to 128 channels through two-dimensional convolution (2D Convolution). At the same time, the channel dimension is further compressed to 1 through a Maxout layer, and repeated 64 times to obtain X input with 64 channel dimensions. Subsequently, X input is input into a residual module, which is composed of the following steps, which effectively extract local and global features while maintaining the hierarchical consistency of the features through step-by-step standardization, convolution, and attention modules.

[0138] First, the input is standardized, and the input feature X input is standardized to unify the feature scale:

[0139]

[0140] wherein Norm is an instance normalization operation, and then a two-dimensional convolution is used to extract preliminary features:

[0141]

[0142] where Conv2D is a two-dimensional convolution operation, in order to capture multi-scale feature information, three different sizes of convolution kernels are used to further extract features, and three groups of features are spliced to fuse different scale information:

[0143]

[0144] X multi = Contact(X1, X2, X3) (14)

[0145] where X1, X2, X3 are the results after extracting features by different convolution kernels, Contact is a concatenation operation, and then X multi Refer to formula (9), (10) for standardization and two-dimensional convolution compression to get Then use the channel attention mechanism (Efficient Channel Attention, ECA) and the spatial attention mechanism (Spatial Attention, SA) to extract features:

[0146]

[0147] z = σ(Conv1D(y)) (16)

[0148]

[0149] f sa = Conv2D([AvgPool(X ECA ); MaxPool(X ECA )]) (18)

[0150] X SA = X ECA ⊙ σ(f sa ) (19)

[0151] Equations (15) to (17) are the processing flow of the channel attention mechanism: GAP is the global maximum pooling, which averages the spatial dimensions of each channel to generate a global feature representation; Conv1D is a one-dimensional convolution operation; σ is the Sigmoid activation function; ⊙ represents element-wise multiplication. Equations (18) to (19) are the spatial attention mechanism: AvgPool and MaxPool are average pooling and maximum pooling operations, respectively; [·; ·] represents the splicing of the pooling results; X SA is the result after extracting features by spatial attention SA.

[0152]

[0153] Formula (20) adds the output of the attention mechanism to the input features through a residual connection, enhances gradient flow, and avoids the problem of gradient disappearance.

[0154] Finally, the output layer of the model, which adds X output to the input X output_norm , is standardized according to formula (9), and then output is performed through two layers of KAN convolution. KAN convolution (Convolutional Kolmogorov-Arnold Network, ConvKAN) adjusts the weights of the convolution kernel adaptively, so that the network can dynamically adjust its learning process according to the characteristics of the input data, which has a significant advantage in capturing the heterogeneity and complex contact patterns between protein chains.

[0155] X KAN = ConvKAN(X output_norm ) (21)

[0156] ConvKAN is a KAN convolution, and X KAN is the result after KAN convolution processing. The output of formula (21) is processed using a KAN convolution layer and a Sigmoid classifier, and a probability matrix (contact or not) is output. The matrix dimension is L A × L B , which represents the contact probability P contact of each pair of residues:

[0157] P contact = Sigmoid(ConvKAN(X KAN )) (22)

[0158] Further, during the overall training of the model, after calculating the loss, the backpropagation algorithm is used for optimization, the AdamW optimizer is selected for optimization, and the cosine annealing learning rate scheduling is used to dynamically adjust the learning rate.

[0159] (1) Loss function: FocalLoss, used to handle the imbalance of label data:

[0160] Focal Loss = -α(1-p y ′) γ log p y ′ (23)

[0161] Where α is the loss weighting for different class samples, and p y ' is the probability of the model predicting the contact category. γ is an adjustable parameter that controls the degree of weight decay of easy-to-classify samples (not contact). In the present application, α is 0.75 and γ is 2.

[0162] (2) Topnand Top-L / k prediction accuracy:

[0163]

[0164] Top n is the top n contacts with the highest probability value in the predicted inter-chain residue contact probability matrix, and the Top n prediction accuracy is the ratio of the actual contact in the top n contacts with the highest predicted probability value. In the present application, n is 1, 10, 25 and 50. Top-L / k is the top L / k contacts with the highest probability value in the predicted inter-chain residue contact probability matrix, and the Top-L / k prediction accuracy is the ratio of the actual contact in the top L (L is the sequence length of the monomer with the shorter sequence length in the two monomers of the heterodimer) contacts with the highest predicted probability value. In the present application, k is 1, 5 and 10.

[0165] (3) Weighted evaluation parameter W score :

[0166]

[0167] wherein α and β are the weight parameters of the Top n prediction accuracy and the Top-L / k prediction accuracy, respectively, and in the present application, α and β are both set to 1.

[0168] Further, in order to prove the effectiveness of the method proposed in the present application, comparative experiments were conducted on the test data set in Table 1, and the experimental results are shown in Table 2.

[0169] Table 2 Results of the present application compared with the prior art advanced method on the test data set

[0170]

[0171] Table 2 shows the performance of multiple models (DeepComplex, GLINTER, CDPred, DeepInter and AttCON-Hetero) on the test data set in Table 1, and the prediction results evaluate their performance on different indicators, including the prediction accuracy of "top1", "top10", "top25", "top50", "topL / 10", "topL / 5", "topL" and the weighted evaluation index "W score ".

[0172] Based on the table, the following is a comparative analysis, highlighting the performance of AttCON-Hetero (this invention) on the test dataset: Top-1 Prediction Accuracy: AttCON-Hetero's Top-1 prediction accuracy is 30.91%, significantly surpassing other models. CDPred follows with an accuracy of 27.27%. AttCON-Hetero also outperforms DeepInter by 12.73 percentage points in Top-1 prediction accuracy. Top-10 Prediction Accuracy: AttCON-Hetero's Top-10 prediction accuracy is 26.36%, also higher than all other models. CDPred ranks second with an accuracy of 21.09%. AttCON-Hetero also outperforms DeepInter by 8.36 percentage points in Top-10 prediction accuracy. Top-25 Prediction Accuracy: AttCON-Hetero's Top-25 prediction accuracy is 23.05%, also higher than all other models. CDPred ranks second with an accuracy of 18.55%. Meanwhile, AttCON-Hetero outperformed DeepInter by 5.96 percentage points in Top 25 prediction accuracy. In Top-50 prediction accuracy, AttCON-Hetero achieved 20.55, again ranking first among all models. This demonstrates that AttCON-Hetero maintains high accuracy even when expanding its prediction range, effectively ranking within a larger prediction set.

[0173] Top-L / 10, Top-L / 5, and Top-L prediction accuracy: AttCON-Hetero achieved prediction accuracy of 24.80%, 21.61%, and 15.66% at these scales, respectively, outperforming CDPred by 5.49, 3.61, and 1.93 percentage points. These results further demonstrate that AttCON-Hetero performs excellently across all accuracy ranking metrics, both within smaller and larger prediction ranges.

[0174] Weighted evaluation parameter W score AttCON-Hetero's weighted evaluation parameter W score AttCON-Hetero scored 162.94, significantly higher than CDPred's 135.53 and DeepInter's 119.51. This means that AttCON-Hetero demonstrates stronger overall capabilities and superiority in the comprehensive evaluation indicators.

[0175] Figure 5The display is the predicted heterodimer interchain residue contact probability map of the sample with PDB number 7EW0 under different methods on the test set, the leftmost is the real contact map, where the yellow part is the contact, and the right three pictures are the predicted contact probability maps of DeepInter, CDPred and AttCON-Hetero respectively, the lighter the color, the greater the contact probability, it can be seen from the figure that AttCON-Hetero is closest to the real contact map, thereby having the highest prediction accuracy. Table 2, Figure 5 It is fully proved that the application has strong prediction ability on the heterodimer interchain residue contact.

[0176] The specific embodiments of the application are described in detail above in combination with the drawings, but the application is not limited to the above embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the purpose of the application.

Claims

1. A method for predicting interchain residue contacts in heterodimers, characterized in that, The method comprises the following steps: S1, screening heterodimer structure information Heterodimer-PDB from a protein database; S2, the data obtained in step S1 is divided into training set, validation set and test set, according to the heterodimer structure information, the Heterodimer-PDB structure data is preprocessed, the structure information MonomerA-PDB / MonomerB-PDB and sequence information MonomerA-FASTA / MonomerB-FASTA of the two monomers of the heterodimer are extracted, and the sequence length is L A and L B , and directly spliced into a joint sequence FASTA with a length of L A +L B ; according to the heterodimer structure information Heterodimer-PDB, the Euclidean distance of the amino acids between the two monomers of the heterodimer is calculated, and the contact matrix between the heterodimers is calculated. S3, using the two monomer sequence information MonomerA-FATSA / MonomerB-FASTA generated in step S2, searching the Uniclust30 sequence database using the HHblits multiple sequence alignment search algorithm to obtain multiple sequence alignment information MSAs of the two monomers, and matching the MSAs of the two monomers using TaxID species classification numbers, and finally generating a joint multiple sequence alignment Paired MSAs; S4, using the joint multiple sequence alignment information Paired MSAs generated in step S3, using two protein language models ESM-MSA and ESM2 to extract row attention features respectively; using the CCMpred algorithm to extract the coevolution information features in the Paired MSAs; using the structure information of the two monomers of the heterodimer in step S2 to calculate the distance map within the amino acid chain of the monomer respectively; using the joint sequence FASTA obtained in step S2, using the PSI-BLAST tool to generate a position-specific scoring matrix (PSSM); S5, fusing the various types of features generated in step S4 to generate a fused feature matrix, and improving the diversity of the input features; S6, inputting the features processed in step S5 into a model AttCON-Hetero composed of a deep residual network ResNet and a KAN (Kolmogorov-Arnold Networks) convolutional network, wherein the model uses an efficient channel attention ECA and a spatial attention SA module to capture local and global dependent features of the feature matrix, and uses a KAN convolution with self-adaptive convolution kernel weight to output a predicted contact probability matrix, and finally outputs a probability matrix of amino acid contact between the two chains of the heterodimer; S7, based on the network output of step S6, using the contact matrix Inter-ChainContact between the two monomers of the heterodimer in step S2 as a label for training, then using the validation set in step S2 to verify the performance of the model obtained by each batch training, and saving the model parameters with the best indicators obtained on the validation set as the weight of the final classification model; Finally, verify the performance of the saved model on the test set in step S2.

2. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, S1 screened the PDB database for heterodimers that met the following criteria: (1) the aggregation state should be heterodimer; (2) the biological assembly should contain only two single chains; (3) the resolution of the X-ray structure should be better than (4) the length of the single chains should be greater than 50 and the total length of the complex should be less than 500; (4) the heterodimers were clustered using a 40% sequence similarity threshold and the best resolved heterodimer from each cluster was selected.

3. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, S2 extracts monomer structure information MonomerA-PDB / MonomerB-PDB and sequence information of the heterodimer, and calculates the contact matrix between the two monomers of the heterodimer, comprising the following steps: S2.1, extract monomer information and sequence information of the heterodimer: Extracting structure information and sequence information of two chains of heterodimers: For qualified dimers, extract the three-dimensional structure information of the two chains from their PDB files, save them as Monomer-PDB A and Monomer-PDB B respectively, and extract their amino acid sequence information MonomerA-FATSA and MonomerB-FASTA, and combine MonomerA-FATSA and MonomerB-FASTA into a sequence FASTA; During the extraction process, remove the residues with amino acid label UNK; S2.2, calculating the Euclidean distance and contact matrix between the amino acid residues of the two chains of the heterodimer: For residues i and j, i from monomer A and j from monomer B, two types of distance calculation rules are defined: i∈MonomerA,j∈MonomerB (1) In formula (1), d is the shortest heavy-atom distance, where r a and r b are the three-dimensional coordinates of the heavy atoms in residues i and j, respectively, and ||r a - r b || is the Euclidean distance of the shortest heavy atom between residue i and residue j. In formula (1), d is the shortest heavy-atom distance, where r a and r b are the three-dimensional coordinates of the heavy atoms in residues i and j, respectively, and ||r a - r b || is the Euclidean distance of the shortest heavy atom between residue i and residue j. In formula (1), d In formula (2) is the distance between the beta carbon atom in residue i and residue j (C β -distance), is the C β atom coordinate of residue i, or is less than or equal to 1, otherwise 0, and saved as Inter-Chain Contact (HeavyAtom) and Inter-Chain Contact (C β ), respectively. The matrix size is L A x L B , which is used as the label of the model for comparison of experimental results.

4. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, S3 specifically includes the following steps: S3.1, using HHblits to search the Uniclust30 database to obtain MSAs: Using the monomer sequence information generated in step S2.1, use the HH-suite 3.30 tool to preliminarily search the Uniclust30 database, and generate MSA information for the two monomers respectively; S3.2, matching the MSAs of the two monomers to generate a joint multiple sequence alignment (Paired MSAs): Match the MSAs of monomer A and monomer B using the species classification number TaxID of each monomer, ensuring that the two monomers come from the same biological species or related homologous species; According to the TaxID information, integrate the matching results of monomer A and monomer B to generate a joint multiple sequence alignment Paired MSAs; This joint MSA contains homologous sequence information from the two monomers, providing high-quality input data for subsequent feature extraction.

5. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, S4 The detailed steps for extracting multiple features using the joint multiple sequence alignment Paired MSAs and other feature sources are as follows: S4.1, extract row attention feature F ESM-MSA and F ESM2 F ESM-MSA i.e. ESM-MSA attention features: Based on the joint multiple sequence alignment information (Paired MSAs) generated in step S3, use the ESM-MSA model to extract the row attention features in the sequence to generate the row attention feature matrix ESM-MSARowAttention; F ESM2 ESM2 row attention feature: using joint sequence FASTA, extracting row attention features through the ESM2 model, generating a row attention feature matrix ESM2 RowAttention; S4.2, extract co-evolution information feature F CCMpred , F DCA-APC and F DCA-DI F CCMpred That is, CCMpred coevolution features: based on Paired MSAs, the probability distribution of coevolution between sequences is calculated using the CCMpred algorithm to generate a diagonally symmetric coevolution matrix; F DCA-APC and F DCA-DI : further based on the Direct Coupling Analysis (DCA) method, using the CCMpred algorithm to generate a DCA feature matrix DCA-APC applying Average Product Correction (APC) and a DCA feature matrix DCA-DI without correction, respectively, which reveal the strength and pattern of sequence co-evolution, for describing the interaction between residues; F DCA-DI = DCA (4) DCA in formula (3) is the direct coupling analysis result, Row Mean, Col Mean and Overall Mean are the row mean, column mean and average of all elements of the DCA matrix respectively; S4.3, extracting two monomer chain internal distance features F intra Based on the two monomer structure information MonomerA-PDB / MonomerB-PDB generated in step S2, calculate the Euclidean distance matrix Intra-Chain Distance between the amino acids within the two monomers; For each pair of residues (a, b) within the monomer chain, its distance can be represented as: In Equation (5), (x a ,y a ,z a ) and (x b ,y b ,z b ) are the three-dimensional coordinates of the a-th and b-th residues, respectively; the intra-chain distance matrix generated by monomer A is D A , and the intra-chain distance matrix generated by monomer B is D B , which can be represented as follows: S4.4, generating a position-specific scoring matrix F PSSM The position-specific scoring matrix (PSSM) of the sequence is calculated by the PSI-BLAST tool using the joint sequence FASTA in step S2; the PSSM captures the substitution patterns and conservation information of amino acids through multiple iterations of searching the UniProt database, and the generated matrix can be used for feature enhancement and model optimization; f in equation (7) ij is the frequency of amino acid j at position i, b j is the background frequency.

6. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, The detailed steps of S5 for generating fusion features are as follows: The seven features generated in step S4 are concatenated in the channel dimension to generate the fusion feature F contact ; F contact = [F ESM-MSA , F ESM2 , F CCMpred , F DCA-APC , F DCA-DI , F intra , F PSSM ] (8) F in equation (8) contact The dimension of F is (L A +L B )×(L A +L B )×C, where L A , L B are the sequence lengths of monomer A and monomer B, respectively, and C represents the number of channels after fusion.

7. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, S6 Network model construction based on deep learning and protein language model includes the following steps: The feature matrix F generated in step S5 is input into an AttCON-Hetero model based on a deep residual network (ResNet) and KAN (Kolmogorov-Arnold Networks) convolution for feature learning, and the dimension of the feature matrix is (L contact +L A )×(L B +L A )×C. B Instance normalization is used to standardize the input feature matrix to reduce the distribution deviation of the input features; two-dimensional convolution is used to compress the C feature channels; meanwhile, the Maxout layer is used to further extract the channel features to obtain X input ; then X input is input into a residual module, each residual module is composed of the following steps, and local and global features are effectively extracted through step-by-step standardization, convolution and attention modules while maintaining the consistency of the feature hierarchy; First, the input features X input are standardized to unify the feature scale: Where Norm is the instance normalization operation, and then a two-dimensional convolution is used to extract preliminary features: Where Conv2D is a two-dimensional convolution operation, in order to capture multi-scale feature information, three different size convolution kernels are used for further extraction, and the three groups of features are spliced to fuse different scale information: X multi = Contact(X1, X2, X3) (14) Wherein X1, X2, X3 are the results after different convolution kernels extract features, Contact is a concatenation operation, and then X multi Refer to formulas (9), (10) to continue the standardization process and two-dimensional convolution compression Then use the channel attention mechanism ECA and the spatial attention mechanism SA to extract features: z = σ(Conv1D(y)) (16) f sa = Conv2D([AvgPool(X ECA ). MaxPool(X ECA )]) (18) X SA = X ECA σ(f sa ) (19) Equations (15) to (17) are the processing flow of the channel attention mechanism: GAP is the global maximum pooling, which averages the spatial dimensions of each channel to generate a global feature representation; Conv1D is a one-dimensional convolution operation; σ is the Sigmoid activation function; and ⊙ represents element-wise multiplication; X ECA is the result after the channel attention ECA extracts the features; and equations (18) to (19) are the spatial attention mechanism: AvgPool and MaxPool are the average pooling and maximum pooling operations, respectively; and [·; ·] represents the concatenation of the pooling results; X SA is the result after the spatial attention SA extracts the features. Formula (20) adds the output of the attention mechanism to the input features through a residual connection to enhance gradient flow and avoid the problem of gradient disappearance; Finally is the output layer of the model, will X output Refer to formula (9) for standardization, get X output_norm , then output through two layers of KAN convolution layer, KAN convolution ConvKAN adjusts the weight of the convolution kernel adaptively, so that the network can dynamically adjust its learning process according to the characteristics of the input data; X KAN = ConvKAN(X output_norm ) (21) ConvKAN is a KAN convolution, X KAN is the result after KAN convolution processing; using a KAN convolution layer and a Sigmoid classifier on the output of formula (21), output a probability matrix (contact or not), the output matrix dimension is L A x L B , represents the contact probability P contact of each pair of residues. P contact = Sigmoid(ConvKAN(X KAN )) (22).

8. The method of predicting interchain residue contacts of a heterodimer according to claim 1, wherein, S7 The steps of model training process, model performance evaluation and model testing are as follows: S7.1 The use of loss function in model training and the optimization of optimizer The contact prediction matrix output in step S6 is taken as the output of the network, and the heterodimer inter-chain residue contact matrix obtained in step S2 is used as the label label for model training; The training process is optimized by the FocalLoss loss function, using the AdamW optimizer combined with the cosine annealing learning rate schedule; S7.2 Performance evaluation and model selection After each batch of training, the model is evaluated for performance using the validation set provided in step S2, by computing Topn accuracy, Top-L / k accuracy and weighted evaluation metric W score to evaluate the model performance, and the model parameters that perform best on the validation set are ultimately selected.

Images