Graph neural network interpretability method based on key sub-graph retrieval

Abstract

The invention discloses a graph neural network interpretability method based on key sub-graph retrieval, and provides a brand-new global and non-model-dependent key sub-graph retrieval method for serving decision interpretation aiming at the problem that an existing graph neural network (GNNs) interpretability strategy often neglects the importance of sub-graphs. According to the strategy, based on node representation generated by GNNs, a core sub-graph influencing graph label attributes can be accurately identified. And through prediction experiments on a BA3 data set and a Mutageicity data set and visual comparison evaluation with other mainstream interpretation strategies, the superior accuracy and excellent generalization of the method provided by the invention are verified.

Description

Technical Field

[0001] The present invention relates to the field of drug design and neural network intelligent recognition, and in particular to a graph neural network interpretability method based on key subgraph retrieval. Background Art

[0002] Graph Neural Networks (GNNs) are widely used in many fields such as social network analysis, physical system simulation and prediction, traffic network optimization, and natural language processing (NLP) due to their excellent performance in complex relationship capture, representation learning, generalization ability, and multi-task processing. And with the joint research of scholars in machine learning and the medical industry, it is found that GNNs are very good at tasks such as graph classification, node classification, and link prediction, and show a high adaptability to molecular property prediction in pharmacology.

[0003] Traditional drug development methods generally rely on tedious experiments and trial and error processes, which are not only inefficient but also full of risks. However, the rise of GNNs technology has opened up a new path for the prediction of drug molecular properties. The molecular properties of drugs, such as toxicity and efficacy, are actually rooted in the exquisite combination of the structure of molecular active functional groups. In the process of drug design and development, GNNs can play its absolute advantages in graph classification, node classification, and link prediction, and show extremely high adaptability to molecular property prediction research in pharmacology.

[0004] Although the practice of graph neural networks (GNNs) in medicine has achieved remarkable results, its "black box" characteristics have caused academic circles to question the rationality and legitimacy of its decision-making, which has greatly affected the credibility of the model in key areas and restricted its further development. [7] With the widespread use of machine learning technology, people from all walks of life are paying more and more attention to understanding the logic behind model decisions, rather than simply pursuing prediction accuracy. The interpretability of graph neural networks has gradually attracted attention.

[0005] The interpretability of GNNs is to make the internal operation mechanism of the network transparent, and then clarify the mapping relationship between the model when processing input and output. Interpretability plays a vital role in establishing system trust, ensuring controllability, discovering and improving model errors, avoiding ethical issues, and meeting regulatory requirements. Especially in the field of drug design, the interpretability of the model reveals the subtle connection between the key components of the drug and its mechanism of action, which is of great significance for improving the accuracy of drug design, improving screening efficiency, and ensuring the reliability of drug efficacy.

[0006] At present, the mainstream algorithms of GNNs' interpretability strategies are mainly classified into the following four categories:

[0007] First, gradient / feature-based interpretation strategies: represented by SA, GuidedBP, CAM, and Grad-CAM, these strategies show the diversity of these strategies in implementation details, especially the customized processing of the back-propagation mechanism and feature fusion technology, by differentially designing the back-propagation process and hidden layer feature integration strategy. However, the interpretation of this type of interpretation strategy is local, and as the data dimension increases, its computational cost increases sharply. Among them, feature-based interpretation methods require human intervention when encountering certain unexplainable features, and there is a possibility of introducing subjective errors.

[0008] Second, perturbation-based explanation strategies: They aim to explore the specific impact of small changes in input data on model output in applications of deep graph models. For example, GNNExplainer12 significantly improves the accuracy of explanations by maximizing the mutual information between GNNs predictions and the possible subgraph structure distribution. PGExplainer parses the prediction mechanism by learning approximate discrete masks of edges, providing a global perspective explanation for trained GNNs, achieving a dual improvement in understanding breadth and consistency. ZORRO uses a greedy algorithm to select nodes with maximized scores during the iterative process of node identification, and replaces the remaining node values ​​with random noise to evaluate the consistency between the new prediction and the initial prediction. GraphMask is trained by minimizing the divergence index of prediction differences over the entire dataset, thereby obtaining a classifier that can determine whether removing a specific edge will affect the predicted output. In addition to having the same problems as gradient / feature-based explanation strategies, this type of strategy is extremely dependent on the choice of perturbations. The range of perturbations is usually limited, which may cause the strategy to fail to capture the behavior of the model in extreme or abnormal situations, and small perturbations may lead to significant changes in the prediction results, which is contrary to the expectations of the model in practical applications.

[0009] Third, agent-based explanation strategy: The core of this strategy is to simplify the complex nonlinear relationships in the graph model and provide instance-level explanations. Specifically, GraphLime introduces the LIME algorithm into the graph field and deeply analyzes the role of node features in the classification process. It uses HSIC Lasso to accurately fit the model, thereby providing an intuitive explanation for GNN predictions based on feature weights. However, this strategy focuses too much on node features and relatively ignores the importance of connections between nodes in the prediction process. In contrast, RelEx adopts a different strategy. It combines agent models with perturbation techniques, simulates local scenes by sampling to generate connected subgraphs, and uses perturbation techniques to parse the logic behind the prediction. This strategy can consider the interactions between nodes more comprehensively, but the limitations of RelEx in context understanding and remote relationship extraction restrict its effectiveness and credibility in practical applications. This limitation is particularly evident when dealing with unstructured, complex or highly specialized text data. CXPlain is a causal model that does not depend on the structure of the original model and is suitable for black box models, but it has a huge demand for data and high computational costs, making it less friendly to non-professionals. Another strategy worth noting is PGM-Explainer, which forms a rich local knowledge base by changing node features and tracking the impact of these changes on prediction results. This strategy is highly general, but there is still room for improvement in considering the value of side information.

[0010] Fourth, the decomposition-based graph model interpretation strategy: This strategy aims to construct an appropriate score distribution criterion to ensure that the decomposed items can be accurately reconstructed into the original prediction value. Such strategies include Layerwise Relevance Propagation (LRP), Excitation BP and GNN-LRP. These strategies have their own characteristics in terms of operating logic and interpretation goals: LRP focuses on calculating node importance indicators to enhance the credibility of the interpretation; Excitation BP emphasizes providing explanations by calculating neuron probability equivalence; and GNN-LRP deeply explores the graph walk process and is closely integrated with GNNs. However, in practical applications, the decomposition-based graph model interpretation strategy may lead to interpretation errors due to over-simplification of the model, and the computational complexity is high and the versatility is insufficient.

[0011] like Figure 1As shown in the figure, in the process of using the GNNs model for prediction, it is the key substructure that can really determine the prediction result. Even if GNNs obtains the correct result through irrelevant information and achieves the expected performance, the credibility of its model is questioned. The above-mentioned various GNNs interpretation strategies are basically focused on the characteristic analysis of nodes and edges, ignoring the core role of key substructures in the decision-making process. In the task of predicting the properties of drug molecules, the core substructure of the molecule often carries the key information of the molecular activity. Therefore, when optimizing the GNN interpretation strategy, identifying and analyzing these key subgraphs can achieve a deep understanding of the model prediction logic. At the same time, most of the existing interpretation methods are related to the model. When facing new tasks or data sets, customized adjustments and optimizations are required. This highly dependent setting limits the versatility of existing strategies. Therefore, it is necessary to develop more flexible and more generalizable interpretation strategies to adapt to the ever-changing scientific research and application environment. Summary of the invention

[0012] In view of the problems existing in the above-mentioned prior art, the present invention provides a graph neural network interpretability method based on key subgraph retrieval, which can accurately obtain the keyword graph that affects the label attributes, improve the interpretability of GNNs in the task of drug molecule property prediction, and the strategy is independent of the model, does not rely on a specific internal structure or data set, and does not require repeated GNNs training, further improving the interpretability and credibility of GNNs in the field of drug design.

[0013] The method comprises the following steps:

[0014] Step 1: Use the BA3 dataset and the Mutagenicity dataset to train the graph neural network model respectively;

[0015] Step 2: For the drug molecule graph g to be designed, all graphs matching its labels are screened out in the BA3 dataset and the Mutagenicity dataset after prediction by the graph neural network model trained in step 1, and these graphs are used as candidate graph sets for graph g;

[0016] Step 3: Use the trained graph neural network model as a node representation extractor, and input the graph g and its candidate graph set into the trained graph neural network model to generate node representations;

[0017] Step 4: Retrieve key subgraphs from the candidate graph set based on the K-nearest neighbor algorithm and output the key subgraphs;

[0018] Step 5: Compare the output key subgraph with the drug molecule graph for interpretable explanation.

[0019] Furthermore, the specific steps of step 1 include:

[0020] Step 1.1: Preprocess the BA3 dataset and divide the preprocessed dataset into training set, test set and validation set;

[0021] Step 1.2: Preprocess the Mutagenicity dataset and divide the preprocessed dataset into training set, test set and validation set;

[0022] Step 1.3: Use the BA3 and Mutagenicity training sets to train the graph neural network models respectively. During the training process, Adam is used as the optimizer, and the learning rate is dynamically adjusted with the learning rate decay strategy. The dropout regularization method is used to process the nodes.

[0023] Step 1.4: Test the trained graph neural network model using the data in the BA3 and Mutagenicity test sets;

[0024] Step 1.5: Verify the performance of the tested graph neural network model on the validation sets of BA3 and Mutagenicity.

[0025] Furthermore, the specific steps of step 1.1 include:

[0026] Step 1.1.1: Construct node features based on the node index and role information of the graph;

[0027] Step 1.1.2: Randomly perturb the node features;

[0028] Step 1.1.3: Set the edge feature to a vector of all 1s;

[0029] Step 1.1.4: Convert the graph labels into PyTorch tensors and adjust the dimensions.

[0030] Step 1.1.5: Extract the node location information from the original data and save it;

[0031] Step 1.1.6: Construct the torch_geometric data object, which contains node features, edge features, edge indexes, edge labels, graph labels, and true masks;

[0032] Step 1.1.7: Divide the BA3 dataset into training, validation, and test sets in a ratio of 2:1:1 and save them in PyTorch’s data format.

[0033] Furthermore, the specific operation steps of step 1.2 include:

[0034] Step 1.2.1: Read the edge index information from the Mutagenicity_A.txt file, construct edge features by mapping chemical valences to integer values, and perform one-hot encoding on edge labels, convert them into corresponding feature representations, and store them as PyTorch tensors;

[0035] Step 1.2.2: Read the node label information from the Mutagenicity_node_labels.txt file, construct node features by mapping the labels of chemical elements to integer values, and perform one-hot encoding on the node labels, convert them into corresponding feature representations, and store them as PyTorch tensors;

[0036] Step 1.2.3: Read the graph label information from the Mutagenicity_graph_labels.txt file and load it as a one-dimensional array using numpy. Convert the graph label information to the PyTorch tensor format and ensure that the data type is long.

[0037] Step 1.2.4: Construct a graph data object, which includes node features, edge features, edge indexes, edge labels, and graph labels;

[0038] Step 1.2.5: Read the graph indicator information from the Mutagenicity_graph_indicator.txt file, divide the dataset into different subgraphs according to the graph indicator information, randomly shuffle the subgraph list, and divide the dataset into training set, validation set, and test set in a ratio of 2:1:1.

[0039] Furthermore, the specific operation steps of step 4 include:

[0040] Step 4.1: Obtain the node representation of graph g and the node representation of the candidate graph;

[0041] Step 4.2: Take graph g as the test graph, use Euclidean distance to measure the similarity between any two points in the test graph and the candidate graph, and obtain the similarity matrix S between graph g and each candidate graph;

[0042] Step 4.3: Based on the similarity matrix S, select the top K key subgraphs g of graph g through a voting mechanism sub .

[0043] Furthermore, the specific steps of step 4.2 include:

[0044] Step 4.2.1: Calculate the similarity s between each node in graph g and each node in all candidate graphs in turn by the following formula: ij :

[0045]

[0046] Among them, is the node representation of the i-th node for testing Figure X 1, and is the node representation of the j-th node of the X2 candidate graph;

[0047] Step 4.2.2: Combine all the similarities s calculated between graph g and each candidate graph ij to form a similarity matrix S.

[0048] Furthermore, the specific steps of Step 4.3 include:

[0049] Step 4.3.1: According to the similarity s in each similarity matrix S ij , select the K nodes with the smallest similarity s ij , and denote these K nodes as set K;

[0050] Step 4.3.2: Compare the similarity s of each node in set K ij with the similarity threshold s_threshold respectively. If s ij < s_threshold, then select this node, and set both the i-th row and the j-th column of the matrix S to which this node belongs to infinity, and enter Step 4.3.3;

[0051] Step 4.3.3: Vote on the selected nodes, and calculate the voting weight of each node using the following formula:

[0052]

[0053] where both α and β are hyperparameters; Con i represents the number of connected components of node i;

[0054] and Con i = ∑ j∈K A ij

[0055]

[0056] where A ij represents the adjacency matrix of node i;

[0057] Step 4.3.4: Output the top K nodes with the highest voting weight as the key subgraph.

[0058] Therefore, by adopting the above method for interpretable graph neural network based on key subgraph retrieval, the present invention has the following beneficial effects:

[0059] First, the present invention improves the accuracy and robustness of interpretation through global perspective retrieval. It considers the overall data distribution and the relationship between candidates, rather than just local data points or individual candidates, by voting from all candidates and giving explanatory conclusions based on the voting results. This global idea can better capture the characteristics and structure of the overall data, better deal with noise and uncertainty in the data, and make the final explanatory results more comprehensive and reliable.

[0060] Second, the present invention uses the powerful representation ability of GNN to consider the similarity of graphs. It converts the candidate graph into node representations by using the trained GNN model, and uses these representations to calculate the similarity with the target graph. Since GNN can effectively capture the semantic information of nodes in the graph structure, the similarity calculation based on GNN representation can more comprehensively compare the structure and features of the graph, thereby improving the accuracy and reliability of the similarity.

[0061] Third, in response to the problem that current explanation strategies often focus on isolated points or edges to describe GNN predictions, ignoring the importance of subgraph connectivity, resulting in inconsistency with facts, the present invention proposes a key subgraph retrieval method that considers the connectivity between substructures. In addition to using node distance to measure subgraph similarity, it also takes into account the connectivity of nodes, which can accurately capture key structures that affect predictions and improve the ability to identify functional groups in molecules.

[0062] Fourthly, the present invention also sets a similarity threshold so that only candidate images with a distance below the threshold are selected for interpretation, thereby avoiding interference of irrelevant candidate images on the interpretation result and improving the accuracy and reliability of the interpretation.

[0063] In summary, the present invention explains the prediction results of the model in an intuitive way by retrieving the candidate graph that is most similar to the graph to be explained and accurately identifying the key subgraphs therein, thereby enhancing the understanding of the decision focus and process of the model. Moreover, the explanation strategy proposed in the present invention does not depend on a specific model and can be flexibly applied to various GNN models. The superiority of this method, as well as its excellent interpretability and generalization, are demonstrated through experimental results, thus providing a new perspective and method for enhancing the interpretability of GNNs in the task of drug molecule property prediction.

[0064] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 Internal decision logic for the GNNs prediction model.

[0066] Figure 2 Flowchart retrieved for key subgraphs.

[0067] Figure 3 This is a diagram of the GNN model training process and technical framework.

[0068] Figure 4 This is the Acc curve of the BA3 train data set.

[0069] Figure 5 This is the loss curve of the BA3 train data set.

[0070] Figure 6 This is the Acc curve of the BA3 validation dataset.

[0071] Figure 7 This is the loss curve of the BA3 validation data set.

[0072] Figure 8 This is the Acc curve of the Mutagenicity train data set.

[0073] Fig. 9 This is the loss curve of the Mutagenicity train data set.

[0074] Fig.10 This is the Acc curve of the Mutagenicity validation dataset.

[0075] Fig.11 This is the loss curve of the Mutagenicity validation dataset.

[0076] Fig.12 First visualization of the mutagenicity explanatory subgraph.

[0077] Fig.13 Second visualization of the Mutagenicity explanatory subgraph.

[0078] Fig.14 Comparison of Mutagenicity explanatory subgraph visualization results. DETAILED DESCRIPTION

[0079] In the description of the present invention, it is also necessary to explain that, unless otherwise clearly specified and limited, these embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention. In addition, it should be understood that after reading the content taught by the present invention, those skilled in the art make various changes or modifications to the present invention, and these equivalent forms also fall within the scope limited by the claims attached to the application.

[0080] In the field of bioinformatics, molecular structures are cleverly mapped into graph models, where atoms become nodes and chemical bonds become edges. Each node atom carries characteristics such as type, charge, and mass; the edges represent different types of chemical bonds and their characteristics, highly summarizing the microscopic world of complex biological molecules. In drug prediction tasks, the core substructure in a molecule often determines the specific properties of the chemical functional group, such as toxicity, activity, or solubility.

[0081] See attached Figure 1-14 Shown

[0082] The present invention provides a graph neural network interpretability method based on key subgraph retrieval, and its workflow is

[0083] like Figure 2 As shown in the figure, it can be seen that it includes the following steps:

[0084] Step 1: Use the BA3 dataset and the Mutagenicity dataset to train the GNN model respectively;

[0085] Step 2: For the drug molecule graph g to be designed, after prediction by the GNN model trained in step 1, all graphs matching its labels are screened out in the BA3 dataset and the Mutagenicity dataset, and these graphs are used as candidate graph sets for graph g;

[0086] Step 3: Input the graph g and its candidate graph set into the trained graph neural network model, and use the trained graph neural network model as a node representation extractor to generate node representations;

[0087] Step 4: Retrieve key subgraphs from the candidate graph set based on the K-nearest neighbor algorithm and output the key subgraphs;

[0088] Step 5: Compare the output key subgraph with the drug molecule graph for interpretable explanation.

[0089] According to the above steps, each part of the method is described in detail below.

[0090] 1. Graph Neural Network Model (GNN) Training

[0091] In the training phase, the present invention selected the BA3 (BA3-Motif) data set and the Mutagenicity data set as model inputs to pre-train the GNN model. Among them, the BA3-Motif database is based on the Barabasi-Albert (BA) graph, and is additionally expanded through any of the three structures of house, cycle, and grid. The formed pattern or sequence is repeated in the biological sequence, which is of great significance to biological functional studies such as pattern discovery, function prediction, and gene regulation. The Mutagenicity data set is mainly used in toxicology research, public health risk assessment, drug development, and biomarker discovery. In the Mutagenicity data set, each molecule is based on a carbon ring structure, among which molecules containing NOz and NI atomic groups are judged to be mutagenic.

[0092] (1) BA3 dataset training

[0093] In the training of the BA3 dataset, the node embedding layer of GNN uses the Linear function to nest a linear transformation to convert the node from the original 4-dimensional features to a 64-dimensional feature vector, that is, self.node_emb=Limnear(4,64). This is to map the original node features to a higher-dimensional space to meet the input requirements of the model and capture more complex feature relationships.

[0094] The convolution layer uses LEConv (Local Extrema Convolution) convolution layer, which is a graph convolution network that can capture local extrema. LEConv is used to calculate the cluster fitness score φ of each cluster i :

[0095]

[0096] Where N(i) represents G c The domain of the i-th node in i A represents the average intra-domain distance in the field, W1, W2, W3 represent learnable parameters, and σ(·) represents the activation function. ij Represents the adjacency matrix of nodes i and j:

[0097]

[0098] The number of input and output channels of each convolution layer of LEConv is 64. A ReLU activation function is added after each convolution layer to increase the nonlinearity of the model.

[0099] The global pooling layer combines the feature vectors of all nodes into a graph-level feature vector through global average pooling (global_mean_pool). This is achieved by simply averaging the features of all nodes, so that a feature vector that represents the entire graph can be obtained.

[0100] The prediction layer uses two fully connected layers and a Softmax function to convert the graph-level feature vector into a prediction result. The two fully connected layers are used to further process the graph-level features, and the Softmax function is used to output the prediction probability of each category. The Sofmax formula is as follows:

[0101]

[0102] Among them, Z i is the output value of the i-th node, and C is the number of output nodes, that is, the total number of categories.

[0103] The Softmax function can convert the output values ​​of multiple categories into a probability distribution between [0, 1] and 1.

[0104] (2) Mutagenicity dataset training

[0105] A different training strategy is used for the Mutagenicity dataset. In the node embedding layer, the original 14-dimensional features of the nodes are converted into 32-dimensional feature vectors, i.e. self.node_emb = Linear(14,32). In addition, since the edges of the Mutagenicity dataset also contain 3-dimensional features, the edge features are also embedded in the same way to meet the data input requirements of the convolutional layer.

[0106] In the convolution layer, the graph representation learning network GINEConv is used, which is a convolutional structure based on the graph isomorphism network. The graph isomorphism network GIN is the most classic graph representation learning network at present. Its principle is as follows:

[0107] Let X (the set of finite multisets) be countable, then there exists a function f: X→R n For any real number x∈X, c∈X, we have:

[0108] h(c,x)=(1+ε)·f(c)+∑ x∈X f(x) (4)

[0109] Among them, h(c,x) is unique for each combination of (c,x), that is, an injective aggregation function. At the same time, any function g can be decomposed into:

[0110] g(c,x)=φ(1+ε)·f(c)+∑ x∈Xf(x) (5)

[0111] Thus, we get the aggregation method of GIN, and then introduce the multi-layer perceptron MLP to learn the functions φ and f:

[0112]

[0113] The global pooling layer and prediction layer are processed in the same way as the BA3 dataset.

[0114] In the training process of both datasets, Adam optimizer is used. Adam’s first-order moment vector m t , second-order moment vector v t The calculation formula is as follows:

[0115] m t =β1m t-1 +(1-β1)g t (7)

[0116] v t =β2v t-1 +(1-β2)·diag(g t 2 ) (8)

[0117] in,

[0118] m0=0 and v0=0;

[0119] β1 represents the exponential decay rate of the first-order moment estimate (usually set to 0.9);

[0120] β2 represents the exponential decay rate of the second-order moment estimate (usually set to 0.999);

[0121] g t represents the gradient at time step t,

[0122]

[0123] At the initial stage of iteration, m t and v t There is a shift towards the initial value. Therefore, a bias correction is made to the first and second order momentum:

[0124]

[0125] This is used to update the Adam optimizer parameters:

[0126]

[0127] Among them, η represents the learning rate. As a key hyperparameter, it determines the size of the step taken towards the minimum value of the loss function; ε represents the smoothing term, which is used to avoid the denominator being zero and is usually taken as 1e-8.

[0128] In addition, the learning rate decay strategy (Learning RateScheduler) is used to dynamically adjust the learning rate η. After each training cycle, the performance of the model on the validation set is evaluated in time, and the model parameters are adjusted based on the evaluation results. The mode is selected as min. The hyperparameters of this decay strategy include factor, patience, and min_lr. Factor is the decay factor used to adjust the learning rate. If the performance indicator on the validation set does not improve within consecutive patience epochs, it will be decayed using the decay factor factor, but the adjustment of η will not be lower than the lower limit of the learning rate min_lr.

[0129] 2. Generation of GNN node representation

[0130] In a graph neural network (GNN), we usually need to generate a vector representation for each node, which will contain information about the node. During the model training process, the node representation is dynamically adjusted. As each batch of data is passed, the node representation is constantly adjusted and updated based on the current features and information of neighboring nodes. This dynamic adjustment enables the model to better capture the structural information and feature distribution in the graph data, thereby improving the generalization ability and performance of the model. As the training progresses, the learning process of the node representation will eventually converge to a stable state after multiple rounds of iterations, and the model performance will also stabilize and no longer change significantly.

[0131] For the BA3 dataset, when generating node representations, the input node features are first embedded through a linear layer to convert the feature dimensions into the node representation dimensions (64 dimensions) required by the model, and dropout regularization is used to avoid overfitting in the generation process. That is, during the training process, dropout randomly sets the outputs of some nodes in the neural network to zero, and discards these units with a certain probability, so that only some nodes participate in the propagation of information each time it is updated, thereby reducing the existence dependency between units and the complexity of the model.

[0132] Secondly, use multi-layer LEConv to aggregate and update node representation layer by layer. Each layer of LEConv will propagate and aggregate information based on neighbor node information and edge attributes, and combine ReLU activation function for nonlinear transformation. Through the above steps, node representation can integrate the information of surrounding nodes and capture the local structural information of graph attributes. In the process of information aggregation, in order to prevent overfitting, dropout regularization is applied to node representation again, and finally a reliable node representation is generated.

[0133] For the Mutagenicity dataset, its node representation generation process is similar to that of the BA3 dataset, except that the convolutional layer used in BA3 is LEConv, which uses the local extreme value information of node features to update the node representation. The convolutional layer used in Mutagenicity is GINEConv, which is a convolutional layer based on graph isomorphism. It captures the structural and local information of the graph by considering the features of neighboring nodes and incorporating them into the node representation.

[0134] 3. Key subgraph retrieval based on K-nearest neighbor algorithm

[0135] The overall idea is to use the trained GNN to generate node representations for the graph for the label prediction task of a certain data set, and use the similarity of the node representations of the graphs with the same prediction results to get the explanation results. The following are the specific steps:

[0136] (1) Generation of candidate graph sets

[0137] In the process of retrieving key subgraphs, it is particularly important to build a comprehensive and accurate set of candidate graphs. If the data set has n categories, then n sets of candidate graphs need to be built, each set corresponding to one category. Then select candidate graphs from the graphs with correct prediction results to ensure the accuracy and reliability of the candidate graphs. This means that the candidate graphs are formed by graphs that are successfully predicted by the model, so the features and structural layouts in the candidate graphs are more likely to be similar to the features and structures of the real graphs, and these candidate graphs have higher interpretability. The candidate graph collection helps to decode the decision logic of complex models and uncover the deep laws and mysteries hidden behind the data.

[0138] (2) Similarity measurement of node representation

[0139] In the process of key subgraph retrieval, it is necessary to first calculate the similarity between the candidate graph and the target graph in order to find the most relevant candidate graph and the substructure with the highest similarity. The trained GNN is used as the extractor of node representation. The present invention uses the K-nearest neighbor algorithm (KNN) to find key substructures.

[0140] For graphs g1 and g2, the trained GNN is used as a node extractor to obtain the node representations X1 and X2 corresponding to the two graphs. The dimension of X1 is (M, D), where M is the number of nodes in graph g1 and D is the feature dimension of each node. The dimension of X2 is (N, D), where N is the number of nodes in graph g2. The two vectors capture the feature information of each node, such as the attributes of the node, the information of neighbor nodes, etc.

[0141] The similarity between node i of graph g1 and node j of graph g2 is measured by Euclidean distance, which is calculated as follows:

[0142]

[0143] in, represents the node representation of the i-th node in graph g1 containing k features, Represents the node representation of the j-th node in graph g2 containing k features.

[0144] The Euclidean distance calculation results of each node in graph g1 and graph g2 are stored in the matrix S. The dimension of matrix S is (M, N). Each value s in the matrix ij Represents the similarity between node i of graph g1 and node j of graph g2. The smaller the value, the more similar they are, and vice versa.

[0145] (3) Key subgraph retrieval strategy

[0146] For the test graph g, firstly, all the graphs that match the label of graph g and are predicted correctly by GNN are selected in the training set and validation set, and these graphs are used as candidate graphs of graph g. Then, graph g is compared with all the candidate graphs. Figure 1 And as the input of the trained GNN model, the node representation of the entire graph is input. During the training process, GNN learns the representation of nodes by considering the neighbor information of nodes and the topological structure of the graph. This node representation method captures the semantic information of nodes in the graph structure and can better interpret the relationship between nodes.

[0147] According to formula (12), the Euclidean distance between each node in graph g and each node in each candidate graph is calculated. Graph g and each candidate graph will calculate an S matrix, which stores the Euclidean distance s between each node in graph g and each node in the candidate graph. ij Since the points with shorter distances are more similar, the values ​​of the corresponding positions on the S matrix are smaller. Therefore, the most similar substructure between graph g and the candidate graph can be determined based on the S matrix, and this substructure is the key factor that leads to the same labels of the two graphs, that is, the key substructure.

[0148] Since there are many candidate graphs, not every candidate graph is useful. Therefore, in order to select objects that are beneficial to the results from a large number of candidate graphs, the present invention defines the following voting mechanism for candidate graphs:

[0149] ① For the graph g in the test set, calculate the S matrix of the test graph g and each candidate graph one by one, and select the K nodes with the smallest distance from all matrices, that is, the first K s ij The smallest node. For the selected s ij , then the i-th row and j-th column of the matrix S are set to infinity to prevent the node i of the graph g from being selected repeatedly.

[0150] ② Vote for the selected K nodes. ij The smaller it is, the higher the similarity between points i and j is, and the greater the possibility of it being used as a guide for GNN prediction results, so the weight obtained should be higher. The weight calculation formula for node i in the definition graph g is:

[0151]

[0152] Among them, α represents a hyperparameter;

[0153] ③ Since the core element that affects the properties of molecules is the functional group, which has a specific connectivity structure, and the current strategy focuses on isolated points or edges to predict GNN, ignoring the importance of subgraph connectivity, resulting in inconsistency with the facts. Therefore, the present invention incorporates connectivity into the voting weight of the node. First, the adjacency matrix of graph g is calculated to store the connectivity status of graph g. Assuming that graph g contains M nodes, the dimension of the adjacency matrix A is (M, M), and its adjacency matrix is ​​shown in formula (2).

[0154] After selecting the K nodes with the smallest distance in the candidate graph, the selected nodes are recorded as set K, and then the number of connections between each node in set K and the other K-1 nodes is calculated. The number of connections between node i in set K is represented by Con i To indicate that Con i It can be calculated by the adjacency matrix, and its calculation formula is:

[0155] Con i =∑ j∈K A ij (14)

[0156] ④ Considering the distance between nodes and the connectivity of nodes at the same time, the weight calculation formula shown in formula (13) is updated to:

[0157]

[0158] Among them, α and β are hyperparameters.

[0159] ⑤ Since the number of candidate graphs may be extremely large, even if two graphs have the same label, but if their similarities are very low, it indicates that the correlation between them is not significant. In this case, continuing to use such candidate graphs for interpretation may lead to incorrect results. Therefore, to determine the value of K, this strategy introduces a similarity threshold. Only when s ij < s_threshold, the candidate graph is used for interpretation, so as to screen out the candidate graphs with higher correlation with the current graph, thereby improving the accuracy and reliability of the interpretation.

[0160] 4. Hyperparameter Tuning

[0161] As can be seen from the above, there are three hyperparameters α, β, and s_threshold in the retrieval strategy proposed by the present invention. Therefore, the present invention uses the Bayesian optimization algorithm to tune the hyperparameters of the model. By adjusting the hyperparameters, we can make the model better fit the data, reduce overfitting, and improve the robustness of the model. Bayesian optimization for parameter tuning is a parameter optimization method based on the Gaussian process and Bayes' theorem. Bayesian optimization is carried out by means of surrogate optimization, that is, by using a surrogate function to fit the objective function, continuously sampling in the hyperparameter space, and updating the surrogate function.

[0162] First, define a reasonable search range for the hyperparameters α, β, and s_threshold, so as to effectively reduce the time cost of parameter tuning. Then set the objective function, whose return value is based on the recall metric on the test dataset (recall), and use 1 - recall as the return value of the objective function. The Bayesian optimization algorithm will dynamically adjust the values of the hyperparameters according to the previous evaluation results, and select the combination of hyperparameters that is most likely to lead to performance improvement in each iteration, ensuring that the model can continuously move towards the goal of performance improvement during the optimization process.

[0163] Embodiment

[0164] To verify the effectiveness of the method proposed by the present invention, the following experiments are carried out.

[0165] 1. Experimental Environment

[0166] The experimental environment configuration in this embodiment is shown in Table 1.

[0167] Table 1 Experimental Environment Configuration

[0168] Configuration Information operating system Ubuntu 20.04.2 Development Language Python 3.9.18 frame PyTorch1.13.1+cu116 CPU Intel Xeon Gold 6248R CPU @ 3.00 GHz GPU GeForceRTX3090(24G) Main Python Libraries torchgeometric2.4.0

[0169] 2. Dataset Preprocessing

[0170] (1) BA3 Dataset Preprocessing

[0171] Load the original graph dataset file BA-3motif.npy and preprocess the original data as follows:

[0172] ①Node feature construction: Based on the node index and role information of the graph, the node features are constructed. These features are integrated into a 4-dimensional vector to accurately describe the different attributes of the node.

[0173] ②Node feature perturbation: The node features are randomly perturbed. A parameter called alpha is introduced to precisely control the degree of perturbation.

[0174] ③ Edge feature setting: In the BA3 dataset, edges mainly carry topological information, which is relatively unimportant in this strategy. To simplify the model input, the edge feature is set to a vector of all 1s.

[0175] ④ Label processing: According to the model input requirements, the graph labels are converted into PyTorch tensor form, and necessary dimension adjustments are made.

[0176] ⑤Position information extraction: extract the node location information from the original data and save it.

[0177] ⑥ Graph data object construction: Construct a torch geometric Data object, which contains key information such as node features, edge features, edge indexes, edge labels, graph labels, and true masks.

[0178] ⑦ Dataset division: The BA3 dataset is divided into training set, validation set, and test set according to the ratio of 2:1:1, and saved in PyTorch data format.

[0179] (2) Mutagenicity dataset preprocessing

[0180] Read the Mutagenicity raw data file, which contains the following key information: edge index, label, graph indicator, label, and node label. The Mutagenicity preprocessing process is as follows:

[0181] ① Edge feature construction: Extract edge index information from the `Mutagenicity_A.txt` file and construct edge features by mapping chemical valence to integer values. Perform one-hot encoding on edge labels, convert them into their corresponding feature representations, and save them in PyTorch tensor format.

[0182] ②Node feature construction: Read the node label information from the `Mutagenicity_node_labels.txt` file, convert it into an integer value using the mapping relationship of chemical elements to form node features. Similarly, perform one-hot encoding on the node labels to obtain the corresponding feature representation and save it as a PyTorch tensor.

[0183] ③ Label processing: Load the graph label information from the `Mutagenicity_graph_labels.txt` file and convert it into a one-dimensional array using numpy. Subsequently, convert the label information into a PyTorch tensor with a data type of 1ong.

[0184] ④ Graph data object construction: Based on the torch_geometric library, construct the Data object to integrate information such as node features, edge features, edge indexes, edge labels, and graph labels.

[0185] ⑤ Dataset partitioning: Read the graph indicator information in the `Mutagenicity_graph indicator.txt` file, and determine the graph to which different nodes belong based on this information. Based on the graph indicator information, the dataset is partitioned into multiple subgraphs, each of which contains specific nodes and corresponding edges. The dataset partitioning and storage format are consistent with BA3.

[0186] 3. Experimental performance indicators

[0187] Given a test set containing a test graph g, this strategy is used to obtain all the key subgraphs g that play a decisive role in the properties of g. sub .

[0188] In order to verify the effectiveness of the proposed method and the decisiveness of the extracted sub-graph g, the extracted sub-graph g sub Re-input into the GNN model and perform the prediction task independently to verify g sub Whether it is possible to reproduce results that are completely consistent with or highly close to the original full-image g prediction.

[0189] (1) Accuracy

[0190] The experiment uses the Accuracy metric to measure the explanatory power of subgraphs, that is, to calculate the proportion of correctly predicted instances based on extracted subgraphs to the total number of test graphs.

[0191] (2) Recall@5

[0192] Accuracy is not enough to fully evaluate the explanatory power of the model. For example, when GNNs capture the features of key subgraphs due to wrong information, their prediction results may not conform to the generally accepted explanations of humans. Therefore, this paper introduces the Recall metric to quantify the model's ability to accurately identify and match key subgraphs in the test set environment.

[0193] Specifically, the Recall@5 algorithm calculates the probability that the model successfully retrieves the true label subgraph nodes under the condition of limiting the subgraph scale (that is, the number of nodes in all predicted graphs does not exceed 5). Recall@5 focuses on the more microscopic, easy-to-understand and explain subgraph structure, effectively eliminating the noise that may be introduced by large-scale complex subgraphs, ensuring the effectiveness of the model while also ensuring its rationality.

[0194] 4. Loss function setting

[0195] For the BA3 dataset and the Mutagenicity dataset, the present invention uses CrossEntropyLoss as the loss factor. The CrossEntropyLoss triggering condition of BA3 is: once the verification loss does not show a significant downward trend within 10 consecutive epochs, the system will decay the learning rate to 80% of the current value, while ensuring that the learning rate will not be lower than 1e-4 to achieve a balance between the learning rate and the risk of overfitting; in the training process of the Mutagenicity dataset, a phased training strategy is adopted. In the first 10% of the epochs, all the data of the training set divided based on the Mutagenicity dataset are directly globally optimized, and then the dynamic node selection mechanism in the prior art is introduced. At the beginning of each epoch, node screening is performed on each graph in the training set, and only some representative nodes are retained to form subgraphs, and training is carried out on these selected subgraphs. The CrossEntropyLoss triggering condition of Mutagenicity is: when the loss of the verification set does not show a significant decrease for 20 consecutive epochs, the learning rate is reduced by 80% of the current value.

[0196] 5. Experimental results

[0197] (1) GNN training results and analysis

[0198] After extensive hyperparameter tuning, the initial learning rate for both BA3 and Mutagenicity datasets was finally determined to be 1e-3. The prediction performance of BA3 and Mutagenicity datasets on the test set is shown in Table 2:

[0199] Table 2 Prediction performance of BA3 and Mutagenicity on the test set

[0200] Dataset BA3 Mutagenicity Accuracy 0.9925 0.8260 Recall 0.9947 0.8832

[0201] On the BA3 test set, the model achieved an accuracy of up to 99.25%. This result proves that the trained GNN model not only has excellent performance on the training set, but also has strong generalization ability, is capable of extracting node representations, and provides strong support for subsequent key subgraph retrieval.

[0202] TensorBoardx is used to monitor the training process (in the following visualization results generated by TensorBoardx, the light color is the real curve, and the dark color is the curve after setting the smoothness, where Smoothness = 0.6) The performance of BA3 in the training set and the validation set is as follows Figures 4 to 7 As shown in the figure, the accuracy of the training set and the validation set both increase steadily with the increase in the number of iterations. Although the accuracy of the validation set is slightly lower than that of the training set, it just reveals that the model has effectively avoided the risk of overfitting while learning general rules. As the number of iterations accumulates, the loss values ​​of the training set and the validation set both show a rapid decline in the early stage and then tend to be stable, which shows the excellent generalization ability of the GNN model on untouched datasets from a visual perspective.

[0203] The performance of the Mutagenicity dataset in the training set and the validation set is as follows: Figures 8 to 11 As shown in the figure. As can be seen from the figure, as the training deepens, the Acc of the Mutagenicity training set shows a stable trend, and the 1oss value also remains at a low level, which clearly shows that the model is gradually approaching the best performance on the training data. Although the accuracy and loss on the validation set fluctuate slightly, they eventually stabilize, which indicates that the model training has entered a stable stage and there is no obvious sign of overfitting. The slight oscillation may be due to the diversity of the validation set samples and the difference in data distribution. In addition, the accuracy of the Mutagenicity dataset on the test set finally settled at 82.6%, which proves that GNN can accurately extract information closely related to molecular properties from molecular graphs.

[0204] (2) Comparative experiment

[0205] In order to evaluate the performance of this strategy, the present invention compares the performance of this strategy with the current mainstream strategies (SA, Grad-CAM, GNN, CXPlain, PGM) on the Mutagenicity and BA3 datasets.

[0206] The control parameter top_ratio is used in the comparative experiment, and its value range is between [0.1, 1.0], with an interval of 0.1. For each top_ratio value, the integral area (ACC-AUC) of the corresponding subgraph accuracy changing with the dynamic change of top_ratio is calculated. ACC-AUC is the average of 10 top_ratios, and the results are shown in Table 3:

[0207] Table 3 Interpretability indicators on BA3 and Mutagenicity datasets

[0208]

[0209] As can be seen from the above table, this strategy significantly surpasses other mainstream strategies in both ACC and Recall@5. Even when limiting the number of subgraph nodes, this strategy can still maintain excellent recall and accuracy.

[0210] (3) Visualization of prediction results

[0211] In order to deeply and intuitively understand the logic of model prediction, the present invention takes the Mutagemicity dataset as an example, and uses the rdMolDraw2D tool in the RDKit library in the field of chemical informatics for further verification. rdMolDraw2D is a drawing tool in the RDKit library, which is used to draw 2D molecular structure diagrams and performs well in molecular graphics rendering. The algorithm designed by the present invention is nested into the molecular object creation process in rdMolDraw2D, and the styles of atoms and bonds are set in advance. The system will identify and extract the molecular substructures that the model pays full attention to during the analysis process, and highlight the microscopic details behind the model's decisions. In the visualization process, the node and edge information in the molecular structure that the model considers to be crucial to the prediction results is accurately captured, and the complete structural map of the molecule is reconstructed using RDKit, and the key substructures are displayed in a highlighted manner. Fig.12 and Fig.13 In the figure, the highlighted area accurately captures the NH2 atom group, showing the decision logic of the model when analyzing complex molecular structures. These diagrams intuitively reveal the key structural elements that affect the model's judgment, which can help researchers and field experts to more deeply understand the model's decision logic and enhance their trust in the model's performance.

[0212] (4) Comparative evaluation of visualization results of multiple strategies

[0213] The comparative experiment selected a series of representative samples from the test set and compared them with a variety of mainstream interpretation strategies. To ensure the consistency of the comparison, the present invention normalized the scores of all interpretation strategies: under the premise of unifying s_threshold, the key subgraphs of each interpretation strategy were screened out, such as Fig.14 As shown in the figure, this strategy can accurately identify the amino (NH2) functional group that has a decisive influence on the molecular label, while the substructures annotated by other strategies have errors. Compared with the scattered and isolated edges or nodes frequently identified by other strategies, this strategy shows higher structural connectivity when identifying key subgraphs. This is the effective benefit obtained by formula (15), which consolidates the ability of this strategy to accurately capture the key structures that affect prediction, especially in the identification of functional groups in molecules.

[0214] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A graph neural network interpretability method based on key subgraph retrieval, characterized in that: The following steps are involved: Step 1: Use the BA3 dataset and the Mutagenicity dataset to train the graph neural network model respectively; Step 2: For the drug molecule graph g to be designed, all graphs matching its labels are screened out in the BA3 dataset and the Mutagenicity dataset after prediction by the graph neural network model trained in step 1, and these graphs are used as candidate graph sets for graph g; Step 3: Use the trained graph neural network model as a node representation extractor, and input the graph g and its candidate graph set into the trained graph neural network model to generate node representations; Step 4: Retrieve key subgraphs from the candidate graph set based on the K-nearest neighbor algorithm and output the key subgraphs; Step 5: Compare the output key subgraph with the drug molecule graph for interpretable explanation.

2. According to claim 1, a graph neural network interpretability method based on key subgraph retrieval is characterized in that: The specific steps of step 1 include: Step 1.1: Preprocess the BA3 dataset and divide the preprocessed dataset into training set, test set and validation set; Step 1.2: Preprocess the Mutagenicity dataset and divide the preprocessed dataset into training set, test set and validation set; Step 1.3: Use the BA3 and Mutagenicity training sets to train the graph neural network models respectively. During the training process, Adam is used as the optimizer, and the learning rate is dynamically adjusted with the learning rate decay strategy. The dropout regularization method is used to process the nodes. Step 1.4: Test the trained graph neural network model using the data in the BA3 and Mutagenicity test sets; Step 1.5: Verify the performance of the tested graph neural network model on the validation sets of BA3 and Mutagenicity.

3. A graph neural network interpretability method based on key subgraph retrieval according to claim 2, characterized in that: The specific steps of step 1.1 include: Step 1.1.1: Construct node features based on the node index and role information of the graph; Step 1.1.2: Randomly perturb the node features; Step 1.1.3: Set the edge feature to a vector of all 1s; Step 1.1.4: Convert the graph labels into PyTorch tensors and adjust the dimensions. Step 1.1.5: Extract the node location information from the original data and save it; Step 1.1.6: Construct the torch_geometric data object, which contains node features, edge features, edge indexes, edge labels, graph labels, and true masks; Step 1.1.7: Divide the BA3 dataset into training, validation, and test sets in a ratio of 2:1:1 and save them in PyTorch’s data format.

4. The graph neural network interpretability method based on key subgraph retrieval according to claim 2 is characterized in that: The specific steps of step 1.2 include: Step 1.2.1: Read the edge index information from the Mutagenicity_A.txt file, construct edge features by mapping chemical valences to integer values, and perform one-hot encoding on edge labels, convert them into corresponding feature representations, and store them as PyTorch tensors; Step 1.2.2: Read the node label information from the Mutagenicity_node_labels.txt file, construct node features by mapping the labels of chemical elements to integer values, and perform one-hot encoding on the node labels, convert them into corresponding feature representations, and store them as PyTorch tensors; Step 1.2.3: Read the graph label information from the Mutagenicity_graph_labels.txt file and load it as a one-dimensional array using numpy. Convert the graph label information to the PyTorch tensor format and ensure that the data type is long. Step 1.2.4: Construct a graph data object, which includes node features, edge features, edge indexes, edge labels, and graph labels; Step 1.2.5: Read the graph indicator information from the Mutagenicity_graph_indicator.txt file, divide the dataset into different subgraphs according to the graph indicator information, randomly shuffle the subgraph list, and divide the dataset into training set, validation set, and test set in a ratio of 2:1:

1.

5. The graph neural network interpretability method based on key subgraph retrieval according to claim 4 is characterized in that: The specific steps of step 4 include: Step 4.1: Obtain the node representation of graph g and the node representation of the candidate graph; Step 4.2: Take graph g as the test graph, use Euclidean distance to measure the similarity between any two points in the test graph and the candidate graph, and obtain the similarity matrix S between graph g and each candidate graph; Step 4.3: Based on the similarity matrix S, select the top K key subgraphs g of graph g through a voting mechanism sub .

6. A graph neural network interpretability method based on key subgraph retrieval according to claim 5, characterized in that: The specific steps of step 4.2 include: Step 4.2.1: Calculate the similarity s between each node in graph g and each node in all candidate graphs in turn by the following formula: ij : in, is the node representation of the i-th node of the test graph X1, is the node representation of the jth node of the candidate graph of X2; Step 4.2.2: Compare all similarities s calculated between graph g and each candidate graph ij Form a similarity matrix S.

7. The graph neural network interpretability method based on key subgraph retrieval according to claim 6 is characterized in that: The specific steps of step 4.3 include: Step 4.3.1: According to the similarity s in each similarity matrix S ij , select the similarity s ij The smallest K nodes, and these K nodes are recorded as set K; Step 4.3.2: Compare the similarity s of each node in the set K ij with the similarity threshold s_threshold respectively. If s ij < s_threshold, select the node, set both the i-th row and the j-th column of the matrix S to which the node belongs to infinity, and go to Step 4.3.3; Step 4.3.3: Vote for the selected nodes and calculate the voting weight of each node using the following formula: Among them, α and β are hyperparameters; Con i represents the number of connected nodes of node i; And Con i =∑ j∈K A ij Among them, A ij represents the adjacency matrix of node i; Step 4.3.4: Output the top K nodes with the highest voting weights as the key subgraph.

Images