Membrane protein microenvironment function prediction method fusing heterogeneous graph convolution and physicochemical modulation

Abstract

This invention discloses a method for predicting the function of membrane protein microenvironments by integrating heterogeneous graph convolution and physicochemical modulation, belonging to the field of bioinformatics. The method includes: obtaining the amino acid sequence of the membrane protein to be predicted; inputting it into ESM2 to extract an evolutionary feature matrix; constructing an original physicochemical feature matrix by selecting 8-dimensional amino acid physicochemical indicators based on the AAindex database; concatenating the original matrix with the Z-score normalized feature matrix; and performing linear projection dimensionality reduction to obtain initial node features; taking the normalized physicochemical features as conditional variables; mapping them to scaling and translation factors through a multilayer perceptron; and performing element-wise affine transformation on the structural features to achieve physicochemical feature modulation; performing cross-modal semantic alignment with the GO terminology prototype vector through a cross-attention module; and outputting the function prediction result by combining ancestor propagation and hierarchical constraint loss functions. This invention realizes multi-level graph structural representation of membrane protein microenvironment semantics and dynamic modulation of physicochemical features.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, and in particular to a method for predicting the function of membrane protein microenvironments by integrating heterogeneous graph convolution and physicochemical modulation. Background Technology

[0002] Membrane protein function prediction is an important direction in proteomics research. Gene Ontology (GO) functional annotation provides a standardized description of the molecular function, biological process, and cellular component of membrane proteins. Existing methods for membrane protein function prediction typically use protein language models to encode amino acid sequences, extract evolutionary features, and then input them into graph neural networks or fully connected networks for classification and prediction. However, existing methods suffer from the following technical shortcomings in feature representation and graph structure construction.

[0003] Existing methods for predicting membrane protein function typically construct protein graph structures using single adjacency relationships, connecting adjacent residues or spatially close residues with uniform edge weights. However, membrane proteins have complex spatial structures, and the interaction strengths between different residues vary significantly. Strongly contacting residue pairs and weakly contacting residue pairs play different roles in functional transfer. Existing methods do not stratify the strength of spatial contacts, and the graph structure cannot reflect the true interaction strength hierarchy, leading to the dilution of information about key functional sites during feature aggregation.

[0004] Existing methods for predicting membrane protein function typically input physicochemical features directly into the network as part of the initial node features, or combine them with evolutionary features through simple splicing. However, the microenvironmental function of membrane proteins is closely related to their local physicochemical properties (such as hydrophobicity, electrostatics, and dissociation characteristics). Existing methods have not established a dynamic modulation mechanism of physicochemical features on structural features, and cannot adaptively adjust the expression of structural features according to the local physicochemical environment of residues. The interaction between structural features and physicochemical features has not been fully explored.

[0005] Existing methods for predicting membrane protein function typically use a uniform weight matrix to handle all edge relationships in graph convolutional networks, without designing independent aggregation weights for sequence adjacency, strong spatial contact, and weak spatial contact relationships. The information propagation mechanisms differ among these relationship types: sequence adjacency reflects local sequence context, strong spatial contact reflects stable functional site interactions, and weak spatial contact reflects dynamic microenvironmental associations. The current method's uniform aggregation approach leads to information confusion among different relationship types, limiting the expressive power of graph convolution.

[0006] Existing methods for predicting membrane protein function typically treat each GO term as an independent classification target when outputting GO function predictions, without considering the parent-child inheritance relationship in the GO ontology hierarchy. When the model predicts a sub-term as a positive sample but misses its parent term, the prediction result violates the logical consistency of the GO ontology. Existing methods lack a hierarchical constraint mechanism, resulting in insufficient biological interpretability of the prediction results.

[0007] Existing methods for predicting membrane protein function typically use a fixed value (e.g., 0.5) as a global threshold to distinguish between contact and non-contact proteins when determining the contact graph threshold. However, the contact probability distribution of different membrane proteins exhibits an overall shift, with some membrane proteins having generally higher or lower contact probabilities. A fixed threshold leads to abnormal edge density in the graph structure of some membrane proteins, affecting the feature aggregation effect of subsequent graph convolution. Summary of the Invention

[0008] Therefore, embodiments of the present invention provide a method for predicting the function of membrane protein microenvironment by fusing heterogeneous graph convolution and physicochemical modulation.

[0009] A method for predicting the function of membrane protein microenvironments by fusing heterogeneous graph convolution and physicochemical modulation is provided, including the following steps:

[0010] S1. Obtain the amino acid sequence of the membrane protein to be predicted, and input the amino acid sequence into the pre-trained protein language model ESM2 to extract the dimension. Evolutionary trait matrix of 1280 × 10 ,in The length of the amino acid sequence;

[0011] Based on the AAindex database, 8-dimensional amino acid physicochemical indicators were screened, and the amino acid sequences were converted into dimensional [value missing] using a lookup table mapping. ×8 primitive physicochemical characteristic matrix ;

[0012] S2, regarding the original physicochemical characteristic matrix Perform Z-score normalization along the feature dimension, the first The standardized formula for calculating the physical characteristics is: =( - ) / ,in For the first Mean value of physicochemical characteristics under global amino acid statistics For the first The standardized physicochemical characteristic matrix is ​​obtained by taking the standard deviation of the physicochemical characteristics under global amino acid statistics. ;

[0013] The standardized physicochemical feature matrix With the evolutionary feature matrix By concatenating the feature dimensions, we obtain a dimension of... Joint characterization of ×1288 ;

[0014] For the joint characterization Perform linear projection and layer normalization to reduce the dimensionality to 512, and obtain the initial node features. The calculation formula is =LayerNorm( · + ,in It is a learnable weight matrix with dimensions of 1288×512. It is the bias vector;

[0015] S3. Obtain the contact map matrix of the membrane protein. The contact diagram matrix The contact diagram matrix is ​​obtained by predicting the amino acid sequence using a protein structure prediction model. The dimension is × , Contact diagram matrix The Middle The residue and the first Contact confidence between residues , Residue number;

[0016] Statistical contact diagram matrix All conditions are met | - The contact confidence of residue pairs with |>1 forms a contact confidence distribution;

[0017] Sort the contact confidence distribution in ascending order and calculate the first quartile of the contact confidence distribution. and the third and fourth quartiles ;

[0018] The first quartile Determined as a weak contact threshold The third quartile Determined as a strong contact threshold ;

[0019] S4, Based on the weak contact threshold and the strong contact threshold Three types of heterogeneous relation adjacency matrices are derived, including sequential adjacency matrices. Strong spatial contact matrix Contact matrix with weak space ;

[0020] The sequence adjacency matrix Only for those who meet the requirements | - Residues with |=1 are assigned a value of 1 at their positions;

[0021] The strong space contact matrix To satisfy > And| - The residue pair with |>1 is assigned a value of 1 at its position;

[0022] The weak space contact matrix To satisfy < ≤ And| - The residue pair with |>1 is assigned a value of 1 at its position;

[0023] S5, based on the initial node features Using the node features and the three types of heterogeneous relation adjacency matrices as the edge structure, a two-layer heterogeneous relation graph convolutional network is constructed for forward propagation;

[0024] For the target node Features of neighbor nodes are aggregated along sequential adjacency, strong spatial contact, and weak spatial contact relationships, respectively, with each of the three types of relationships employing an independent learnable weight matrix. , , ;

[0025] Aggregation by target node The number of neighboring nodes under the corresponding relationship is the normalization coefficient, which is superimposed on the target node. The self-features are weighted by the self-weight matrix The transformed result is then processed by a nonlinear activation function to obtain the output features of this layer.

[0026] The first layer outputs the first structural feature matrix. The second layer outputs the second structural feature matrix. ;

[0027] S6. Obtain the standardized physicochemical characteristic matrix. The Middle The 8-dimensional physicochemical feature vector corresponding to each residue As a condition variable;

[0028] Will Mapped to a 512-dimensional scaling factor via a first multilayer perceptron The calculation formula is = · + ;

[0029] Will Mapped to a 512-dimensional translation factor via a second multilayer perceptron. The calculation formula is = · + ,in , For learnable weight matrix, , It is the bias vector;

[0030] For the second structural feature matrix The Middle Structural features of each residue Perform an element-wise affine transformation to obtain the modulated features. The calculation formula is = ⊙ + , where ⊙ represents element-wise multiplication;

[0031] S7. Aggregate the modulated features of all residues into a protein characterization, and perform cross-modal semantic alignment between the protein characterization and the GO term prototype vector through a cross-attention module, and calculate the matching score between the membrane protein and each GO term.

[0032] During the training phase, ancestor propagation is performed on the true GO labels of each protein, and the matching score is constrained by a hierarchical constraint loss function to output the membrane protein function prediction results.

[0033] It should be noted that the evolutionary feature matrix in step S1 The extraction method is as follows: the last hidden state of the pre-trained protein language model ESM2 outputs a 1280-dimensional vector for each residue in the amino acid sequence. These 1280-dimensional vectors are then stacked in sequence to obtain a vector with the following dimensions: Evolutionary trait matrix of 1280 × 10 ;

[0034] Contact diagram matrix in step S3 The method for obtaining the amino acid sequence is as follows: the amino acid sequence is input into a protein structure prediction model, the structure prediction model outputs the contact probabilities between all residue pairs, and the contact probabilities are written position by position into a dimension of... × The matrix, Take the first The residue and the first The numerical value of the contact probability between residues;

[0035] Weak contact threshold in step S3 With strong contact threshold The contact confidence distribution of each membrane protein was calculated. The reason for this implementation is that the contact probability distribution of different membrane proteins has an overall shift. When using fixed values ​​to divide strong and weak contacts, the strong spatial contact edges of some membrane proteins are too dense or too sparse, and the graph structure cannot reflect the true interaction intensity level.

[0036] The threshold is set using the quartile of the contact confidence distribution itself. The division of the three types of relationship edges is adaptively adjusted according to the contact confidence distribution of the samples, and is not sensitive to outlier confidence.

[0037] The projection dimension 512 in step S2 is determined as follows: the joint representation dimension 1288 is compressed to 512 by one linear projection. The number of parameters in the linear projection layer is 1288×512+512=659968. The hidden dimension of the heterogeneous graph convolution is consistent with the projection dimension, avoiding the quadratic increase in the number of parameters caused by the high-dimensional representation directly entering the graph convolution.

[0038] Scaling factor in step S6 Perform channel-by-channel scaling and translation factor on each feature channel of the structural features. Baseline shift is performed on structural features, and both use physicochemical features as conditional signals to enable the local hydrophobic environment, electrostatic environment and dissociation environment of the residues to participate in the correction of semantic features;

[0039] In step S7, the ancestor propagation and the hierarchical constraint loss function work together to ensure that the prediction results remain logically consistent along the parent-child hierarchy of the GO ontology.

[0040] Beneficial effects

[0041] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following:

[0042] This invention obtains the contact graph matrix of membrane proteins, statistically analyzes the contact confidence distribution of residue pairs with a residue spacing greater than 1, and calculates the first and third quartiles in ascending order as weak and strong contact thresholds. Based on the weak and strong contact thresholds, it derives three types of heterogeneous relationship adjacency matrices: a sequence adjacency matrix, a strong spatial contact matrix, and a weak spatial contact matrix. This solves the technical problems of existing technologies that use fixed global thresholds to divide contact and non-contact, and do not consider the overall shift in the contact probability distribution of different membrane proteins. It achieves the technical effect of adaptively determining thresholds through contact confidence distribution and constructing a three-type heterogeneous relationship graph structure that reflects the true interaction strength level.

[0043] This invention uses initial node features as node features and three types of heterogeneous relation adjacency matrices as edge structures to construct a two-layer heterogeneous relation graph convolutional network for forward propagation. For the target node, neighbor node features are aggregated along sequential adjacency relations, strong spatial contact relations, and weak spatial contact relations, respectively. Each of the three types of relations uses an independent learnable weight matrix, with the number of neighbor nodes as the normalization coefficient. The results of the feature transformation are superimposed and output through nonlinear activation. This solves the technical problem of using a unified weight matrix to handle all edge relations and the confusion of information of multiple types of relations in the prior art. It achieves the technical effect of realizing the differential feature aggregation of multiple types of relations and improving the expressive power of graph convolution through independent weight matrices.

[0044] This invention takes the physicochemical feature vectors of each residue in the standardized physicochemical feature matrix as condition variables, maps them to scaling and translation factors through a multilayer perceptron, and performs element-wise affine transformation on the structural features output by the second layer. This allows the local hydrophobic environment, electrostatic environment, and dissociation environment of the residue to participate in the correction of semantic features. This solves the technical problems in the prior art where the physicochemical features have not established a dynamic modulation mechanism for structural features and the interaction between the local physicochemical environment and structural features has not been fully explored. It achieves the technical effect of realizing adaptive conditional modulation of structural features by physicochemical features through feature modulation network and enhancing the semantic expression of microenvironment. Attached Figure Description

[0045] Figure 1 A flowchart provided for an embodiment of this application. Detailed Implementation

[0046] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0047] like Figure 1 As shown, the method for predicting membrane protein microenvironment function by integrating heterogeneous graph convolution and physicochemical modulation includes the following steps:

[0048] S1. Obtain the amino acid sequence of the membrane protein to be predicted, and input the amino acid sequence into the pre-trained protein language model ESM2 to extract the dimension. Evolutionary trait matrix of 1280 × 10 ,in The length of the amino acid sequence;

[0049] Based on the AAindex database, 8-dimensional amino acid physicochemical indicators were screened, and the amino acid sequences were converted into dimensional [value missing] using a lookup table mapping. ×8 primitive physicochemical characteristic matrix ;

[0050] S2, regarding the original physicochemical characteristic matrix Perform Z-score normalization along the feature dimension, the first The standardized formula for calculating the physical characteristics is: =( - ) / ,in For the first Mean value of physicochemical characteristics under global amino acid statistics For the first The standardized physicochemical characteristic matrix is ​​obtained by taking the standard deviation of the physicochemical characteristics under global amino acid statistics. ;

[0051] The standardized physicochemical feature matrix With the evolutionary feature matrix By concatenating the feature dimensions, we obtain a dimension of... Joint characterization of ×1288 ;

[0052] For the joint characterization Perform linear projection and layer normalization to reduce the dimensionality to 512, and obtain the initial node features. The calculation formula is =LayerNorm( · + ,in It is a learnable weight matrix with dimensions of 1288×512. It is the bias vector;

[0053] S3. Obtain the contact map matrix of the membrane protein. The contact diagram matrix The contact diagram matrix is ​​obtained by predicting the amino acid sequence using a protein structure prediction model. The dimension is × , Contact diagram matrix The Middle The residue and the first Contact confidence between residues , Residue number;

[0054] Statistical contact diagram matrix All conditions are met | - The contact confidence of residue pairs with |>1 forms a contact confidence distribution;

[0055] Sort the contact confidence distribution in ascending order and calculate the first quartile of the contact confidence distribution. and the third and fourth quartiles ;

[0056] The first quartile Determined as a weak contact threshold The third quartile Determined as a strong contact threshold ;

[0057] S4, Based on the weak contact threshold and the strong contact threshold Three types of heterogeneous relation adjacency matrices are derived, including sequential adjacency matrices. Strong spatial contact matrix Contact matrix with weak space ;

[0058] The sequence adjacency matrix Only for those who meet the requirements | - Residues with |=1 are assigned a value of 1 at their positions;

[0059] The strong space contact matrix To satisfy > And| - The residue pair with |>1 is assigned a value of 1 at its position;

[0060] The weak space contact matrix To satisfy < ≤ And| - The residue pair with |>1 is assigned a value of 1 at its position;

[0061] S5, based on the initial node features Using the node features and the three types of heterogeneous relation adjacency matrices as the edge structure, a two-layer heterogeneous relation graph convolutional network is constructed for forward propagation;

[0062] For the target node Features of neighbor nodes are aggregated along sequential adjacency, strong spatial contact, and weak spatial contact relationships, respectively, with each of the three types of relationships employing an independent learnable weight matrix. , , ;

[0063] Aggregation by target node The number of neighboring nodes under the corresponding relationship is the normalization coefficient, which is superimposed on the target node. The self-features are weighted by the self-weight matrix The transformed result is then processed by a nonlinear activation function to obtain the output features of this layer.

[0064] The first layer outputs the first structural feature matrix. The second layer outputs the second structural feature matrix. ;

[0065] S6. Obtain the standardized physicochemical characteristic matrix. The Middle The 8-dimensional physicochemical feature vector corresponding to each residue As a condition variable;

[0066] Will Mapped to a 512-dimensional scaling factor via a first multilayer perceptron The calculation formula is = · + ;

[0067] Will Mapped to a 512-dimensional translation factor via a second multilayer perceptron. The calculation formula is = · + ,in , For learnable weight matrix, , It is the bias vector;

[0068] For the second structural feature matrix The Middle Structural features of each residue Perform an element-wise affine transformation to obtain the modulated features. The calculation formula is = ⊙ + , where ⊙ represents element-wise multiplication;

[0069] S7. Aggregate the modulated features of all residues into a protein characterization, and perform cross-modal semantic alignment between the protein characterization and the GO term prototype vector through a cross-attention module, and calculate the matching score between the membrane protein and each GO term.

[0070] During the training phase, ancestor propagation is performed on the true GO labels of each protein, and the matching score is constrained by a hierarchical constraint loss function to output the membrane protein function prediction results.

[0071] It should be noted that the evolutionary feature matrix in step S1 The extraction method is as follows: the last hidden state of the pre-trained protein language model ESM2 outputs a 1280-dimensional vector for each residue in the amino acid sequence. These 1280-dimensional vectors are then stacked in sequence to obtain a vector with the following dimensions: Evolutionary trait matrix of 1280 × 10 ;

[0072] Contact diagram matrix in step S3 The method for obtaining the amino acid sequence is as follows: the amino acid sequence is input into a protein structure prediction model, the structure prediction model outputs the contact probabilities between all residue pairs, and the contact probabilities are written position by position into a dimension of... × The matrix, Take the first The residue and the first The numerical value of the contact probability between residues;

[0073] Weak contact threshold in step S3 With strong contact threshold The contact confidence distribution of each membrane protein was calculated. The reason for this implementation is that the contact probability distribution of different membrane proteins has an overall shift. When using fixed values ​​to divide strong and weak contacts, the strong spatial contact edges of some membrane proteins are too dense or too sparse, and the graph structure cannot reflect the true interaction intensity level.

[0074] The threshold is set using the quartile of the contact confidence distribution itself. The division of the three types of relationship edges is adaptively adjusted according to the contact confidence distribution of the samples, and is not sensitive to outlier confidence.

[0075] The projection dimension 512 in step S2 is determined as follows: the joint representation dimension 1288 is compressed to 512 by one linear projection. The number of parameters in the linear projection layer is 1288×512+512=659968. The hidden dimension of the heterogeneous graph convolution is consistent with the projection dimension, avoiding the quadratic increase in the number of parameters caused by the high-dimensional representation directly entering the graph convolution.

[0076] Scaling factor in step S6 Perform channel-by-channel scaling and translation factor on each feature channel of the structural features. Baseline shift is performed on structural features, and both use physicochemical features as conditional signals to enable the local hydrophobic environment, electrostatic environment and dissociation environment of the residues to participate in the correction of semantic features;

[0077] In step S7, the ancestor propagation and the hierarchical constraint loss function work together to ensure that the prediction results remain logically consistent along the parent-child hierarchy of the GO ontology.

[0078] As an optional embodiment, the 8-dimensional amino acid physicochemical indicators in step S1 cover four types of physical constraints;

[0079] The four types of physical constraints include hydrophobicity indicators;

[0080] Indicators of electrostatic field include charge, isoelectric point, and polarity;

[0081] Spatial packing and steric hindrance parameters, including molecular weight and solvent-accessible surface area;

[0082] Local chemical reactivity indicators, including first-order dissociation constant and second-order dissociation constant .

[0083] It should be noted that the screening steps for the 8-dimensional amino acid physicochemical indicators are as follows: entries with missing values ​​are removed from the amino acid index entries included in the AAindex database, and candidate indicators are selected according to the four types of physical constraints.

[0084] Calculate the Pearson correlation coefficient for each pair of candidate indicators, and remove items whose absolute correlation coefficient with the selected indicators is greater than 0.8. 0.8 is the threshold for determining a high degree of correlation in statistics.

[0085] The eight dimensions determined after screening are as follows: hydrophobicity is taken as the Kyte-Doolittle hydrophobicity index, charge is taken as the net charge of the side chain, isoelectric point is taken as the isoelectric point of the amino acid, polarity is taken as the Grantham polarity value, molecular weight is taken as the molecular weight of the residue, and solvent accessible surface area is taken as the statistical value of the average solvent accessible surface area of ​​the residue. Take the dissociation constant of the α-carboxyl group. Take the dissociation constant of the α-amino group;

[0086] The implementation method of the lookup table mapping is as follows: A mapping table of 20 standard amino acids to 8-dimensional values ​​is pre-established. The amino acid sequence is traversed, and the table is looked up position by position. The 8-dimensional values ​​of each position are stacked in sequence order to obtain the dimension. ×8 primitive physicochemical characteristic matrix ;

[0087] The correspondence between the four types of physical constraints and the transmembrane microenvironment of membrane proteins is as follows: hydrophobicity index corresponds to the anchoring stability of transmembrane helices in the lipid core; electrostatic field index corresponds to the electrostatic selectivity of ion channels; spatial stacking and steric hindrance index corresponds to the stacking density of the transmembrane region; and local chemical reactivity index corresponds to the protonation potential of proton pump proteins.

[0088] As an optional embodiment, the mean value in step S2 With the standard deviation It is obtained from global amino acid statistics and remains consistent across the training, validation, and test sets.

[0089] It should be noted that the mean value mentioned above with standard deviation The calculation steps are as follows: Take the first of 20 standard amino acids... Physical index values , This refers to the amino acid type number. Take 1 to 20;

[0090] mean Equal to 20 standard amino acids The arithmetic mean of the dimensional values ​​is calculated as follows:

[0091]

[0092] Standard deviation Equal to 20 standard amino acids Dimensional values ​​relative to the mean The root mean square deviation is calculated as follows:

[0093]

[0094] As an optional embodiment, the first quartile in step S3 With the third quartile The calculation steps are as follows: Let | - The number of residue pairs with |>1 is All contact confidence levels are sorted in ascending order as follows: ≤ ≤…≤ ;

[0095] Ranking of the first quartile =1+( -1)×0.25, ranking of the third and fourth quartiles =1+( -1)×0.75;

[0096] When the rank is not an integer, the first quartile is obtained by performing linear interpolation between the integer part of the rank and the next rank value. and the third quartile ;

[0097] When the rank is an integer, the sort value corresponding to that rank is directly taken as the first quartile. Or the third quartile .

[0098] It should be noted that the calculation steps for the linear interpolation are as follows: Let the ranking be... The integer part is Then the weak contact threshold equal to the first quartile The calculation formula is:

[0099]

[0100] Rank The integer part is Then the strong contact threshold Equal to the third and fourth quartiles The calculation formula is:

[0101]

[0102] The first quartile uses the value at the 25th percentile of all contact confidence levels as the weak contact threshold, and the third quartile uses the value at the 75th percentile as the strong contact threshold. The values ​​of the two thresholds vary with the shape of the contact confidence distribution and are not manually specified.

[0103] As an optional embodiment, in step S4, the sequence adjacency matrix The strong spatial contact matrix and the weak spatial contact matrix All are converted to a sparse coordinate list format for storage;

[0104] The sparse coordinate list format only records the row and column numbers of the residue pairs assigned a value of 1;

[0105] The forward propagation in step S5 performs the collection operation of neighbor node features according to the row number and column number.

[0106] It should be noted that the implementation of the sparse coordinate list format is as follows: for any adjacency matrix, traverse all residue pairs, and assign a value of 1 to the residue index pair. i,j Write the coordinates into the list sequentially;

[0107] Let the sum of the number of edges in the adjacency matrices of the three heterogeneous relation classes be... The storage size of the sparse coordinate list is 2. The storage requirement for a dense adjacency matrix is ​​n integers. × One element;

[0108] During forward propagation, for each edge in the coordinate list i,j Extract the node The features are accumulated to the nodes after being transformed by the weight matrix of the correspondence relationship. The aggregated cache, the computational cost of a single propagation and Proportional;

[0109] Number of edges that satisfy the contact condition With sequence length Approximately linear growth, the storage and computational costs of dense adjacency matrices increase with... The quadratic growth sparse coordinate list format reduces the growth order of storage and computation from quadratic to linear.

[0110] As an optional embodiment, the calculation process for each layer of forward propagation in step S5 is as follows: for the target node With any relation , will node In relationship The features of all neighboring nodes are summed after transformation by the learnable weight matrix corresponding to this relationship, and then divided by the node. In relationship The set of neighboring nodes Number of elements , obtain relationship The normalized aggregation result;

[0111] The normalized aggregation results of sequence adjacency, strong spatial contact, and weak spatial contact are combined with the node... The self-features are weighted by the self-weight matrix The transformed results are summed, and the ReLU activation function is used to output the node features of this layer.

[0112] The heterogeneous relational graph convolutional network has two propagation layers.

[0113] It should be noted that the complete calculation formula for each layer of forward propagation is:

[0114]

[0115] in For the first Layer nodes Features For the first Layer nodes Features It is a set of relations consisting of three types of relations. For nodes In relationship The set of neighboring nodes below, For relationship The corresponding learnable weight matrix, These are self-weight matrices, each with a dimension of 512×512;

[0116] Divide by the number of neighboring nodes The normalization method ensures that the aggregation amplitudes of high-degree nodes and low-degree nodes are consistent, and the characteristic amplitudes of residues in densely contacted regions are not systematically too large.

[0117] The experimental process for determining the number of propagation layers to be two layers is as follows: heterogeneous relational graph convolutional networks with one, two, and three layers are built for comparative training.

[0118] In a single layer, the receptive field only covers first-order neighbors, and long-range spatial interactions cannot reach the target node. The validation set Fmax is lower than that of a two-layer structure.

[0119] In the three-layer structure, the features of each node tend to be consistent after multiple aggregations, and the Fmax of the validation set is also lower than that of the two-layer structure. Therefore, the data of the propagation layer is determined to be two layers.

[0120] The three types of relationships are each assigned an independent weight matrix. The weight matrix corresponding to the sequence adjacency relationship learns the continuous orientation of the peptide backbone, the weight matrix corresponding to the strong spatial contact relationship learns the interaction pattern of high-confidence spatial stacking, and the weight matrix corresponding to the weak spatial contact relationship learns the transmission pattern of long-range weak interactions.

[0121] As an optional embodiment, in step S6, the first structural feature matrix is ​​further processed. Perform the first element-wise affine transformation;

[0122] The first element-wise affine transformation uses independent third and fourth multilayer perceptrons to generate the corresponding scaling and translation factors, respectively.

[0123] The result of the first element-wise affine transformation is used as the input to the second layer of the two-layer heterogeneous graph convolutional network.

[0124] It should be noted that the first, second, third, and fourth multilayer perceptrons are all single-layer linear mappings with an input dimension of 8 and an output dimension of 512. The parameter calculation formula for a single multilayer perceptron is as follows:

[0125]

[0126] The sum of the parameters of the four multilayer perceptrons is 4 × 4608 = 18432;

[0127] The total number of model parameters consists of 659,968 linear projection layers, 2,101,248 heterogeneous graph convolutional layers, 1,050,624 cross-attention modules, and 18,432 multilayer perceptrons, totaling 3,830,272. The proportion of the four multilayer perceptron parameters to the total model parameters is calculated using the following formula:

[0128]

[0129] The first element-wise affine transformation is located between two layers of graph convolution, and the message passing of the second layer is carried out on the features that have been corrected for physicochemical conditions. The local physicochemical environment of the residues participates in the regulation during the neighbor information aggregation stage.

[0130] For the second structural feature matrix The element-wise affine transformation performed is located at the output of the graph convolution, affecting the second structural feature matrix. Perform overall calibration;

[0131] Two modulation locations share a condition variable Each layer uses independent mapping parameters to adapt to the feature distribution of the intermediate and output layers respectively.

[0132] As an optional embodiment, in step S7, the ancestor propagation is implemented as follows: when a protein is labeled as a certain GO term, all ancestor terms of the GO term are traversed along the directed acyclic graph of the GO ontology, and all ancestor terms are marked as positive samples.

[0133] The hierarchical constraint loss function includes hierarchical focus loss and hierarchical consistency penalty term;

[0134] The hierarchical consistency penalty applies a penalty to the portion of a sub-term matching score that is higher than the matching score of its parent term.

[0135] It should be noted that the formula for calculating the hierarchical consistency penalty is:

[0136]

[0137] in t,u In the Go ontology, these are pairs of terms that have a parent-child relationship. As the parent term, For sub-terms, The matching score for the parent term. For the matching score of sub-terms, max(0, - Indicates taking 0 and - The larger value in;

[0138] The hierarchical consistency penalty term takes a non-zero value only when the matching score of the sub-term exceeds the matching score of the parent term, constraining the matching score to be monotonically non-increasing along the GO hierarchy from the root to the leaf.

[0139] Balance weights for each category in hierarchical focus loss The frequency of positive samples by category is calculated using the following formula:

[0140]

[0141] in For category The positive sample frequency in the training set is divided by the sum of the reciprocals of the positive sample frequencies of all classes, and the balance weights are determined by the statistics of the training data.

[0142] The focusing parameter of the hierarchical focus loss was optimally determined to be 2 based on the validation set Fmax during comparative training with values ​​of 1, 2, and 3.

[0143] Ancestor propagation enables the supervision signal of deep sub-terms to be propagated back along the hierarchy to all ancestor terms, and shallow terms are supervised, which together with the hierarchical consistency penalty term constrains the hierarchical logic of the prediction results.

[0144] As an optional implementation, the model is trained using the AdamW optimizer with an initial learning rate of 8e-5;

[0145] The learning rate scheduling employs a combination of a preheating strategy of linear warm-up for the first 5 epochs and a cosine decay strategy.

[0146] Batch size is 24, and mixed precision training is used;

[0147] A weighted attenuation coefficient of 1e-4 is applied to the parameters of the first multilayer perceptron and the second multilayer perceptron. 2. Regularization;

[0148] The early stop mechanism is triggered based on the overall Fmax metric of the validation set.

[0149] It should be noted that the initial learning rate of 8e-5 was determined by learning rate interval testing: the learning rate was increased exponentially from 1e-7 to 1e-2, the training loss was recorded at each level, the curve of training loss as a function of learning rate was plotted, and the learning rate level corresponding to the interval with the largest loss decrease rate was selected. After retesting within this level, it was determined to be 8e-5.

[0150] The batch size of 24 was determined as follows: the memory usage of a single sample was estimated based on the training memory capacity and the maximum sequence length of membrane proteins in the dataset, and the maximum batch size that could be accommodated under memory constraints was taken, which was 24 according to actual measurements.

[0151] The total number of training epochs is 100. The warm-up phase is the first 5 epochs. During the warm-up phase, the learning rate increases linearly with each epoch to the initial learning rate, and then decreases according to a cosine curve.

[0152] The weight decay coefficient 1e-4 was trained on three values, 1e-5, 1e-4, and 1e-3, and then optimally determined according to the validation set Fmax.

[0153] Mixed-precision training uses half-precision floating-point to perform forward and backward propagation, and uses single-precision floating-point to save model parameters;

[0154] The early stopping mechanism is triggered when the overall Fmax of the validation set has no improvement for 10 consecutive epochs.

[0155] As an optional embodiment, in step S7, the GO term prototype vector covers three namespaces: molecular function, biological process, and cellular component.

[0156] The prototype vector of the GO term is obtained by encoding the functional description text of the GO term using a text encoding model;

[0157] The cross-attention module loads the GO term prototype vectors of each namespace, calculates the matching scores, and outputs the prediction results for each namespace.

[0158] It should be noted that the prototype vector of a GO term is constructed as follows: the functional description text of the GO term is taken and encoded into a vector of the same dimension as the protein representation through a text encoding model, which is then used as the prototype vector of the GO term.

[0159] The formula for calculating the cross-attention module is:

[0160]

[0161] Among them, protein characterization is used as a query. Using GO terminology prototype vectors as keys Sum , For the feature dimension, softmax represents the normalization exponential function;

[0162] The GO ontology consists of three independent namespaces: molecular function, biological process, and cellular component. GO term prototype vectors are loaded into each namespace and matching scores are calculated separately. The prediction results of each namespace do not interfere with each other.

[0163] The model evaluation uses four metrics for each namespace: Fmax, AUPR, Smin, and MCC. Fmax is the maximum F1 score obtained by scanning according to the prediction threshold, AUPR is the area under the precision-recall curve, Smin is the semantic minimum distance, and MCC is the Matthews correlation coefficient.

[0164] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A method for predicting the function of membrane protein microenvironment by integrating heterogeneous graph convolution and physicochemical modulation, characterized in that, Includes the following steps: S1. Obtain the amino acid sequence of the membrane protein to be predicted, and input the amino acid sequence into the pre-trained protein language model ESM2 to extract the dimension. Evolutionary trait matrix of 1280 × 10 ,in The length of the amino acid sequence; Based on the AAindex database, 8-dimensional amino acid physicochemical indicators were screened, and the amino acid sequences were converted into dimensional [value missing] using a lookup table mapping. ×8 primitive physicochemical characteristic matrix ; S2, regarding the original physicochemical characteristic matrix Perform Z-score normalization along the feature dimension, the first The standardized formula for calculating the physical characteristics is: =( - ) / ,in For the first Mean value of physicochemical characteristics under global amino acid statistics For the first The standardized physicochemical characteristic matrix is ​​obtained by taking the standard deviation of the physicochemical characteristics under global amino acid statistics. ; The standardized physicochemical feature matrix With the evolutionary feature matrix By concatenating the feature dimensions, we obtain a dimension of... Joint characterization of ×1288 ; For the joint characterization Perform linear projection and layer normalization to reduce the dimensionality to 512, and obtain the initial node features. The calculation formula is =LayerNorm( · + ,in It is a learnable weight matrix with dimensions of 1288×512. It is the bias vector; S3. Obtain the contact map matrix of the membrane protein. The contact diagram matrix The contact diagram matrix is ​​obtained by predicting the amino acid sequence using a protein structure prediction model. The dimension is × , Contact diagram matrix The Middle The residue and the first Contact confidence between residues , Residue number; Statistical contact diagram matrix All conditions are met | - The contact confidence of residue pairs with |>1 forms a contact confidence distribution; Sort the contact confidence distribution in ascending order and calculate the first quartile of the contact confidence distribution. and the third and fourth quartiles ; The first quartile Determined as a weak contact threshold The third quartile Determined as a strong contact threshold ; S4, Based on the weak contact threshold and the strong contact threshold Three types of heterogeneous relation adjacency matrices are derived, including sequential adjacency matrices. Strong spatial contact matrix Contact matrix with weak space ; The sequence adjacency matrix Only for those who meet the requirements | - Residues with |=1 are assigned a value of 1 at their positions; The strong space contact matrix To satisfy > And| - The residue pair with |>1 is assigned a value of 1 at its position; The weak space contact matrix To satisfy < ≤ And| - The residue pair with |>1 is assigned a value of 1 at its position; S5, based on the initial node features Using the node features and the three types of heterogeneous relation adjacency matrices as the edge structure, a two-layer heterogeneous relation graph convolutional network is constructed for forward propagation; For the target node Features of neighbor nodes are aggregated along sequential adjacency, strong spatial contact, and weak spatial contact relationships, respectively, with each of the three types of relationships employing an independent learnable weight matrix. , , ; Aggregation by target node The number of neighboring nodes under the corresponding relationship is the normalization coefficient, which is superimposed on the target node. The self-features are weighted by the self-weight matrix The transformed result is then processed by a nonlinear activation function to obtain the output features of this layer. The first layer outputs the first structural feature matrix. The second layer outputs the second structural feature matrix. ; S6. Obtain the standardized physicochemical characteristic matrix. The Middle The 8-dimensional physicochemical feature vector corresponding to each residue As a condition variable; Will Mapped to a 512-dimensional scaling factor via a first multilayer perceptron The calculation formula is = · + ; Will Mapped to a 512-dimensional translation factor via a second multilayer perceptron. The calculation formula is = · + ,in , For learnable weight matrix, , It is the bias vector; For the second structural feature matrix The Middle Structural features of each residue Perform an element-wise affine transformation to obtain the modulated features. The calculation formula is = ⊙ + , where ⊙ represents element-wise multiplication; S7. Aggregate the modulated features of all residues into a protein characterization, and perform cross-modal semantic alignment between the protein characterization and the GO term prototype vector through a cross-attention module, and calculate the matching score between the membrane protein and each GO term. During the training phase, ancestor propagation is performed on the true GO labels of each protein, and the matching score is constrained by a hierarchical constraint loss function to output the membrane protein function prediction results. Contact diagram matrix in step S3 The method for obtaining the amino acid sequence is as follows: the amino acid sequence is input into a protein structure prediction model, the structure prediction model outputs the contact probabilities between all residue pairs, and the contact probabilities are written position by position into a dimension of... × The matrix, Take the first The residue and the first The numerical value of the contact probability between residues; Weak contact threshold in step S3 With strong contact threshold The contact confidence distribution of each membrane protein was calculated. The reason for this implementation is that the contact probability distribution of different membrane proteins has an overall shift. When using fixed values ​​to divide strong and weak contacts, the strong spatial contact edges of some membrane proteins are too dense or too sparse, and the graph structure cannot reflect the true interaction intensity level. The threshold is set using the quartile of the contact confidence distribution itself. The division of the three types of relationship edges is adaptively adjusted according to the contact confidence distribution of the samples, and is not sensitive to outlier confidence. The projection dimension 512 in step S2 is determined as follows: the joint representation dimension 1288 is compressed to 512 by one linear projection. The number of parameters in the linear projection layer is 1288×512+512=659968. The hidden dimension of the heterogeneous graph convolution is consistent with the projection dimension, avoiding the quadratic increase in the number of parameters caused by the high-dimensional representation directly entering the graph convolution. Scaling factor in step S6 Perform channel-by-channel scaling and translation factor on each feature channel of the structural features. Baseline shift is performed on structural features, and both use physicochemical features as conditional signals to enable the local hydrophobic environment, electrostatic environment and dissociation environment of the residues to participate in the correction of semantic features; In step S7, the ancestor propagation and the hierarchical constraint loss function work together to ensure that the prediction results remain logically consistent along the parent-child hierarchy of the GO ontology.

2. The method for predicting membrane protein function based on heterogeneous relational graph networks and physicochemical modulation according to claim 1, characterized in that, The 8-dimensional amino acid physicochemical indicators in step S1 cover four types of physical constraints; The four types of physical constraints include hydrophobicity indicators; Indicators of electrostatic field include charge, isoelectric point, and polarity; Spatial packing and steric hindrance parameters, including molecular weight and solvent-accessible surface area; Local chemical reactivity indicators, including first-order dissociation constant and second-order dissociation constant .

3. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 2, characterized in that, The mean in step S2 With the standard deviation It is obtained from global amino acid statistics and remains consistent across the training, validation, and test sets.

4. The method for predicting membrane protein function based on heterogeneous relational graph networks and physicochemical modulation according to claim 1, characterized in that, The first quartile in step S3 With the third quartile The calculation steps are as follows: Let | - The number of residue pairs with |>1 is All contact confidence levels are sorted in ascending order as follows: ≤ ≤…≤ ; Ranking of the first quartile =1+( -1)×0.25, ranking of the third and fourth quartiles =1+( -1)×0.75; When the rank is not an integer, the first quartile is obtained by performing linear interpolation between the integer part of the rank and the next rank value. and the third quartile ; When the rank is an integer, the sort value corresponding to that rank is directly taken as the first quartile. Or the third quartile .

5. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 4, characterized in that, In step S4, the sequence adjacency matrix The strong spatial contact matrix and the weak spatial contact matrix All are converted to a sparse coordinate list format for storage; The sparse coordinate list format only records the row and column numbers of the residue pairs assigned a value of 1; The forward propagation in step S5 performs the collection operation of neighbor node features according to the row number and column number.

6. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 1, characterized in that, The calculation process for each layer of forward propagation in step S5 is as follows: for the target node With any relation , will node In relationship The features of all neighboring nodes are summed after transformation by the learnable weight matrix corresponding to this relationship, and then divided by the node. In relationship The set of neighboring nodes Number of elements , obtain relationship The normalized aggregation result; The normalized aggregation results of sequence adjacency, strong spatial contact, and weak spatial contact are combined with the node... The self-features are weighted by the self-weight matrix The transformed results are summed, and the ReLU activation function is used to output the node features of this layer. The heterogeneous relational graph convolutional network has two propagation layers.

7. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 6, characterized in that, In step S6, the first structural feature matrix is ​​also processed. Perform the first element-wise affine transformation; The first element-wise affine transformation uses independent third and fourth multilayer perceptrons to generate the corresponding scaling and translation factors, respectively. The result of the first element-wise affine transformation is used as the input to the second layer of the two-layer heterogeneous graph convolutional network.

8. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 1, characterized in that, In step S7, the ancestor propagation is implemented as follows: when a protein is labeled as a certain GO term, all ancestor terms of the GO term are traversed along the directed acyclic graph of the GO ontology, and all ancestor terms are marked as positive samples. The hierarchical constraint loss function includes hierarchical focus loss and hierarchical consistency penalty term; The hierarchical consistency penalty applies a penalty to the portion of a sub-term matching score that is higher than the matching score of its parent term.

9. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 8, characterized in that, The model was trained using the AdamW optimizer with an initial learning rate of 8e-5. The learning rate scheduling employs a combination of a preheating strategy of linear warm-up for the first 5 epochs and a cosine decay strategy. Batch size is 24, and mixed precision training is used; A weighted attenuation coefficient of 1e-4 is applied to the parameters of the first multilayer perceptron and the second multilayer perceptron.

2. Regularization; The early stop mechanism is triggered based on the overall Fmax metric of the validation set.

10. The method for predicting membrane protein function based on heterogeneous relation graph network and physicochemical modulation according to claim 1, characterized in that, In step S7, the GO term prototype vector covers three namespaces: molecular function, biological process, and cellular component. The prototype vector of the GO term is obtained by encoding the functional description text of the GO term using a text encoding model; The cross-attention module loads the GO term prototype vectors of each namespace, calculates the matching scores, and outputs the prediction results for each namespace.

Images