A method for predicting toxicity of a compound by a graph attention network embedded with molecular isomorphic edge features

Abstract

A method for predicting compound toxicity by embedding molecular isomerization edge features in a graph attention network, comprising the following steps: step 1, preprocessing data to generate a molecular graph; step 2, extracting isomerization edge features of the molecule; step 3, embedding the isomerization edge features into the molecular graph and giving a certain weight; step 4, constructing a graph attention network and modeling using the new molecular graph, and using the trained model for prediction. The present application considers the isomerization edge features in the molecule, models the molecular graph using the graph attention network, and integrates the isomerization edge features into the molecular graph, which better improves the prediction accuracy of the toxicity of compounds such as mutagenicity and cardiotoxicity.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics, specifically relating to a method for predicting compound toxicity using a graph attention network that embeds molecular isomorphic edge features. Background Technology

[0002] Toxicity is often a major cause of failure throughout the entire drug development process, from design to market launch. Today, toxicity evaluation of compounds is considered as important as pharmacodynamic studies. However, traditional drug toxicity testing methods typically require animal models. These animal experiments are not only time-consuming and costly but also raise animal welfare concerns. While high-throughput in vitro screening and other methods have enabled faster evaluation of compound toxicity in recent years, experimental evaluation methods remain costly and time-consuming. A method that predicts compound toxicity using graph attention networks embedded with molecular isomeric edge features, leveraging structural information to predict toxicity, can ensure large-scale toxicity screening before in vitro and in vivo testing. This allows for the early screening of potentially hazardous compounds in drug development, offering a better solution to the problem of toxicity in candidate drugs. Summary of the Invention

[0003] The purpose of this invention is to provide a method for predicting compound toxicity by embedding graph attention networks with molecular isomorphic edge features, which can effectively improve the prediction accuracy of compound mutagenicity, cardiotoxicity and other toxicities.

[0004] To achieve the above objectives, the present invention adopts the following technical solution: a method for predicting compound toxicity through a graph attention network embedding molecular isomorphic edge features, comprising the following steps:

[0005] Step 1: Preprocess the data to generate a molecular diagram;

[0006] Step 2: Extract the isomerism edge features of the molecule;

[0007] Step 3: Embed the heterogeneous edge features into the molecular graph and assign weights;

[0008] Step 4: Construct a graph attention network and model it using the new molecular graph, then use the trained model for prediction.

[0009] Step 1 specifically involves:

[0010] Step 1.1: Preprocess the smiles strings of compounds in the training, testing, and validation sets used for the toxicity prediction model, including desalting, deduplication, removal of mixtures and polymers, removal of compounds containing heavy metals, and removal of compounds that cannot generate molecular diagrams.

[0011] Step 1.2: Initialize the molecular diagram of the compound based on the smiles string.

[0012] Step 2 specifically involves:

[0013] Step 2.1: Starting with a chemical bond in the molecule, find all the triplet pairs in the molecule. A triplet pair includes a neighboring atom, a bond, and a target atom.

[0014] Step 2.2: Determine the molecular isomerism edge based on the elemental composition of the triplet. That is, if there are two or more triplets in which the neighboring atoms and the target atom are the same but the bonds are different, then such bonds constitute the isomerism edge in the molecule.

[0015] Step 3 specifically involves:

[0016] Step 3.1: Generate a heterogeneous edge feature matrix and embed the heterogeneous edge features into the molecular graph;

[0017] Step 3.2: Apply the leaky_relu activation function to the state vectors of the target atom and its neighboring atoms to obtain the heterogeneous edge feature vector and the atom pair feature vector. Then, use the softmax function to assign certain weights to the heterogeneous edge feature vectors. The calculation method is as follows:

[0018] h vu =leaky_relu(W·[g v ,g u (15)

[0019] e vu =leaky_relu(W·[h v ,h u (16)

[0020]

[0021] Among them, g v ,h v It is the state vector of the target atom, g u ,h u is the state vector of the neighboring atoms; W is the trainable weight matrix; h vu It is the feature vector of the heterogeneous edges after activation by the activation function, e vu It is the ordinary edge feature vector after activation by the leaky_relu activation function; N(v) is the set of all neighboring atoms of the target atom; a vu It is the weight of the neighbor atom u relative to the target atom v after being weighted by the softmax function.

[0022] Step 4 specifically involves:

[0023] Step 4.1: Assign the weight a of the neighboring atom u to the target atom v. vu Using the ELU activation function for nonlinear activation (sentence segmentation complete), the context vector C of the target atom v is obtained. v The calculation method is as follows:

[0024]

[0025] Where W is a trainable weight matrix; h u is the state vector of the neighboring atoms; N(v) is the set of all neighboring atoms of the target atom.

[0026] Step 4.2: Neighboring atom u continuously aggregates information for target atom v using the graph attention mechanism. The calculation method is as follows:

[0027]

[0028] Where M k-1 It is a message passing function that iterates k-1 times; It is the state vector of the neighboring atom u after k-1 iterations; N(v) is the state vector of the target atom v after k-1 iterations; N(v) is the set of all neighboring atoms of the target atom. It is the aggregation context feature vector of the target atom after k-1 iterations.

[0029] Step 4.3: Update the state vector of the target atom and obtain the previous state vector of the target atom v. and aggregated context vectors As input, the previous state Update to current state The calculation method is as follows:

[0030]

[0031] GRU k-1 It is an update function that iterates k-1 times; It is the aggregation context feature vector of the target atom after k-1 iterations; It is the state vector of the target atom v after k-1 iterations; It is the state vector of the target atom v after k iterations.

[0032] Step 4.4: Take the state vector of the target atom v after k iterations as input and place it in the output layer of the network to predict the toxicity of the molecule. The calculation method is as follows:

[0033]

[0034] Where W is the dimension reduction matrix set according to the dimension of the atomic feature vector; is the state vector of the target atom v after k iterations; output is the output layer; molFeature is the final prediction result.

[0035] The toxicity prediction model is as follows:

[0036] h vu =leaky_relu(W·[g v ,g u ]) (twenty two)

[0037] e vu =leaky_relu(W·[h v ,h u ]) (twenty three)

[0038]

[0039]

[0040]

[0041]

[0042]

[0043] Among them, the formula for calculating the feature vector of the heterogeneous edge is (8), the formula for calculating the feature vector of the atom pair is (9), the formula for calculating the weight of the neighbor atom to the target atom is (10), the formula for calculating the context feature vector of the target atom is (11), the formula for calculating the aggregation context feature vector of the target atom is (12), the formula for calculating the latest state vector of the target atom is (13), and the formula for predicting the toxicity of the compound based on the latest state vector of the target atom is (14).

[0044] Compared with the prior art, the present invention has the following advantages:

[0045] 1. The modeling process of this invention is simple and clear. It does not require the calculation of redundant molecular fingerprint data. Instead, it uses molecular diagrams with clearer structures for modeling, resulting in higher prediction accuracy for the same compound.

[0046] 2. This invention uses graph attention neural networks to model molecular graphs, which significantly improves the accuracy of model prediction compared to traditional machine learning algorithms.

[0047] 2. In the modeling process, this invention takes into account the isomerism features in the molecule and integrates the isomerism features into the molecular graph, thereby increasing the weight of the isomerism features in the entire molecular graph and improving the accuracy of predicting the mutagenicity, cardiotoxicity and other toxicities of the compound.

[0048] 3. The training set, test set, and validation set data used in the modeling process of this invention come from a wide range of sources, and the model has good generalization ability. Attached Figure Description

[0049] Figure 1 Modeling flowchart;

[0050] Figure 2 ROC plot of compound mutagenicity model on external validation set;

[0051] Figure 3 RMSE plot of the hERG cardiotoxicity model on external validation set 1;

[0052] Figure 4 RMSE plot of the hERG cardiotoxicity model on external validation set 2;

[0053] Figure 5 RMSE plot of the hERG cardiotoxicity model on external validation set 3;

[0054] Figure 6 This is the RMSE plot of the hERG cardiotoxicity model of the present invention on external validation set 4. Detailed Implementation

[0055] The present invention will be further described below with reference to the accompanying drawings and examples.

[0056] Example 1:

[0057] The present invention will be further described below with reference to the accompanying drawings and examples.

[0058] like Figure 1 As shown, the specific embodiment of the present invention is as follows:

[0059] Step 1: For cardiotoxic molecules, 11,512 cardiotoxic molecular compounds were first obtained from the ChEMBL database as the training set, and 1,161 cardiotoxic molecular compounds were selected from the PubChem BioAssay dataset, hERG central database, Open Source Malaria database, and literature as the external validation set. The RDkit chemical toolkit was used for desalting, and the compounds were uploaded to the OCHEM website for deduplication. Compounds that could not generate molecular diagrams were removed, resulting in 9,200 training data and 1,146 external validation data. For mutagenic molecules, mutagenic molecular compounds were obtained from the Ames mutagenicity benchmark dataset as the training set, and mutagenic molecular compounds were obtained from the CCRIS database, the National Toxicology Program NTP database, and the Sunita Institute for Advanced Study in Salmonella Typhi ISSSTY database as the external validation set. Mixtures, polymers, compounds containing heavy metals, and compounds with fewer than 3 carbon atoms were filtered out, and compounds that could not generate molecular diagrams were removed, resulting in 6,506 training data and 1,471 external validation data.

[0060] Step 2: For each smiles string, generate the corresponding molecular diagram. Use the Rdkit chemistry toolkit to find all the bonds in the molecule and the atoms at both ends of the bonds, forming a triplet of neighboring atom, bond, and target atom. Based on all the triplets, determine the isomerism edges in the molecule, that is, triplets where the neighboring atoms and target atoms are the same but the bonds are different. The different bonds constitute the isomerism edges in the molecule.

[0061] Step 3: Initialize the molecular bonds to an all-zero matrix and the heterogeneous edges to 1. Embed the heterogeneous edge features into the molecular graph. Perform nonlinear activation on the state vectors of the target atom and its neighboring atoms using the leaky_relu activation function to obtain the heterogeneous edge feature vector and the atom pair feature vector. Then, use the softmax function to assign certain weights to the heterogeneous edge feature vectors. The calculation method is as follows:

[0062] h vu =leaky_relu(W·[g v ,g u (1)

[0063] e vu =leaky_relu(W·[h v ,h u (2)

[0064]

[0065] Where g v ,h v and g u ,hu These are the state vectors of the target atom and its neighboring atoms, respectively; W is the trainable weight matrix; h vu and e vu These are the heterogeneous edge feature vector and the atom pair feature vector after activation by the leaky_relu activation function, respectively; N(v) is the set of all neighboring atoms of the target atom; a vu The weight is the importance of the neighbor atom u to the target atom v after being weighted by the softmax function.

[0066] Step 4: Assign the importance of neighboring atom u to target atom v, i.e., the weight a. vu The context vector C of the target atom v is obtained by using the ELU activation function for nonlinear activation. v The calculation method is as follows:

[0067]

[0068] Where W is a trainable weight matrix; h u is the state vector of the neighboring atoms; N(v) is the set of all neighboring atoms of the target atom.

[0069] Neighboring atom u continuously aggregates information for target atom v using a graph attention mechanism. The calculation method is as follows:

[0070]

[0071] Where M k-1 It is a message passing function that iterates k-1 times; It is the state vector of the neighboring atom u after k-1 iterations; N(v) is the state vector of the target atom v after k-1 iterations; N(v) is the set of all neighboring atoms of the target atom. It is the aggregation context feature vector of the target atom after k-1 iterations.

[0072] Update the state vector of the target atom and obtain the previous state vector of the target atom v. and aggregated context vectors As input, the previous state Update to current state The calculation method is as follows:

[0073]

[0074] GRU k-1 It is an update function that iterates k-1 times; It is the aggregation context feature vector of the target atom after k-1 iterations; It is the state vector of the target atom v after k-1 iterations; It is the state vector of the target atom v after k iterations.

[0075] The state vector of the target atom v after k iterations is used as input and placed in the output layer of the network to predict the toxicity of the molecule. The calculation method is as follows:

[0076]

[0077] Where W is the dimension reduction matrix set according to the dimension of the atomic feature vector; is the state vector of the target atom v after k iterations; output is the output layer; molFeature is the final prediction result.

[0078] Furthermore, the toxicity prediction model is specifically as follows:

[0079] h vu =leaky_relu(W·[g v ,g u (8)

[0080] e vu =leaky_relu(W·[h v ,h u ]) (9)

[0081]

[0082]

[0083]

[0084]

[0085]

[0086] Formula 8 is used to calculate the feature vector of the heterogeneous edge, Formula 9 is used to calculate the feature vector of the atom pair, Formula 10 is used to calculate the weight of the neighbor atom to the target atom, Formula 11 is used to calculate the context feature vector of the target atom, Formula 12 is used to calculate the aggregation context feature vector of the target atom, Formula 13 is used to calculate the latest state vector of the target atom, and Formula 14 is used to predict the toxicity of the compound based on the latest state vector of the target atom.

[0087] Example 1:

[0088] Given a compound from the CCRIS database, CAS number 820-75-7, to predict its mutagenicity, firstly, a molecular graph is generated based on its smiles string. Then, the Rdkit chemistry toolkit is used to find all isomer edges in the compound. Next, isomer edge features are initialized and embedded into the molecular graph. Based on the molecular graph data, the mutagenicity model of this invention is used for prediction. (Based on the current state) The result of 1 indicates that this compound is mutagenic. The predicted result is consistent with the experimental result.

[0089] Example 2:

[0090] Given a compound from the ChEMBL database, ID: CHEMBL100415, to predict its hERG cardiotoxicity, firstly, a molecular graph is generated based on its smiles string. Then, the Rdkit chemistry toolkit is used to find all isomer edges in the compound. Next, isomer edge features are initialized and embedded into the molecular graph. Based on the molecular graph data, the hERG cardiotoxicity model of this invention is used for prediction. (Based on the current state...) The predicted pIC50 value was 6.89, while the actual pIC50 value was 6.72. The predicted result is close to the experimental result.

[0091] Validation of the effectiveness of this invention: For the prediction of compound mutagenicity, as shown in Table 1, the AUC of this invention on the external validation set is 0.842. This invention was compared with other methods for predicting compound mutagenicity, among which the best AUC on the external validation set in the MutagenPred-GCNNs method is 0.838. In patent CN112466399A, the AUC on the external validation set is 0.797. This invention shows an improvement in the prediction of compound mutagenicity. As shown in Table 2, the RMSE values ​​of this invention on external validation sets 1 to 4 for predicting cardiotoxicity caused by compound hERG were 0.948, 0.908, 0.810, and 1.060, respectively, with a validation set average of 0.9315. This invention was compared with other methods for predicting cardiotoxicity caused by compound hERG. In the TSSF-hERG method, the RMSE values ​​on external validation sets 1 to 4 were 0.876, 0.661, 1.117, and 1.117, respectively, with a validation set average of 0.94275. Although the performance of this invention decreased on validation sets 1 and 2, its performance improved on validation sets 3 and 4. This invention demonstrates improved performance in predicting the cardiotoxicity of compound hERG.

[0092] Table 1 compares the performance of the model with other models on the mutagenicity dataset.

[0093] Methods for predicting the mutagenicity of compounds External validation set AUC This invention 0.842 MutagenPred-GCNNs 0.838 Patent No.: CN112466399A 0.797

[0094] Table 2 compares the performance of the models with other models on the cardiotoxicity dataset.

[0095] Predictive methods for the cardiotoxicity of compound hERG Validation set 1 Validation set 2 Validation set 3 Validation set 4 validation set average This invention 0.948 0.908 0.810 1.060 0.9315 TSSF-hERG 0.876 0.661 1.117 1.117 0.94275 .

Claims

1. A method for predicting compound toxicity using a graph attention network embedding molecular isomerism edge features, comprising the following steps: Step 1: Preprocess the data to generate a molecular diagram; Step 2: Extract the isomerism edge features of the molecule; Step 2 specifically involves: Step 2.1: Starting with a chemical bond in the molecule, find all the triplet pairs in the molecule. A triplet pair includes a neighboring atom, a bond, and a target atom. Step 2.2: Determine the molecular isomerism edge based on the elemental composition of the triplet. That is, if there are two or more triplets in which the neighboring atoms and the target atom are the same but the bonds are different, then such bonds constitute the isomerism edge in the molecule. Step 3: Embed the heterogeneous edge features into the molecular graph and assign weights; Step 3 specifically involves: Step 3.1: Generate a heterogeneous edge feature matrix and embed the heterogeneous edge features into the molecular graph; Step 3.2: Apply the leaky_relu activation function to the state vectors of the target atom and its neighboring atoms to obtain the heterogeneous edge feature vector and the atom pair feature vector. Then, use the softmax function to assign certain weights to the heterogeneous edge feature vectors. The calculation method is as follows: h vu =leaky_relu(W·[g v ,g u ])(1) e vu =leaky_relu(W·[h v h u ])(2) Among them, g v ,h v It is the state vector of the target atom, g u ,h u is the state vector of the neighboring atoms; W is the trainable weight matrix; h vu It is the feature vector of the heterogeneous edges after activation by the activation function, e vu It is the ordinary edge feature vector after activation by the leaky_relu activation function; N(v) is the set of all neighboring atoms of the target atom; a vu It is the weight of the neighbor atom u relative to the target atom v after being weighted by the softmax function; Step 4: Construct a graph attention network and model it using the new molecular graph, then use the trained model for prediction.

2. The method for predicting compound toxicity using a graph attention network embedding molecular isomorphic edge features according to claim 1, wherein step 1 specifically comprises: Step 1.1: Preprocess the smiles strings of compounds in the training, testing, and validation sets used for the toxicity prediction model, including desalting, deduplication, removal of mixtures and polymers, removal of compounds containing heavy metals, and removal of compounds that cannot generate molecular diagrams. Step 1.2: Initialize the molecular diagram of the compound based on the smiles string.

3. The method for predicting compound toxicity using a graph attention network embedding molecular isomorphic edge features according to claim 1, wherein step 4 specifically comprises: Step 4.1: Assign the weight a of the neighboring atom u to the target atom v. vu The context vector C of the target atom v is obtained by using the ELU activation function for nonlinear activation. v The calculation method is as follows: Where W is a trainable weight matrix; h u is the state vector of the neighboring atoms; N(v) is the set of all neighboring atoms of the target atom; Step 4.2: Neighboring atom u continuously aggregates information for target atom v using the graph attention mechanism. The calculation method is as follows: Where M k-1 It is a message passing function that iterates k-1 times; It is the state vector of the neighboring atom u after k-1 iterations; N(v) is the state vector of the target atom v after k-1 iterations; N(v) is the set of all neighboring atoms of the target atom. It is the aggregation context feature vector of the target atom after k-1 iterations; Step 4.3: Update the state vector of the target atom and obtain the previous state vector of the target atom v. and aggregated context vectors As input, the previous state Update to current state The calculation method is as follows: GRU k-1 It is an update function that iterates k-1 times; It is the aggregation context feature vector of the target atom after k-1 iterations; It is the state vector of the target atom v after k-1 iterations; It is the state vector of the target atom v after k iterations; Step 4.4: Take the state vector of the target atom v after k iterations as input and place it in the output layer of the network to predict the toxicity of the molecule. The calculation method is as follows: Where W is the dimension reduction matrix set according to the dimension of the atomic feature vector; is the state vector of the target atom v after k iterations; output is the output layer; molFeature is the final prediction result.

4. The method for predicting compound toxicity using a graph attention network embedding molecular isomorphic edge features according to claim 2, wherein the toxicity prediction model is as follows: h vu =leaky_relu(W·[g v ,g u ])(8) e vu =leaky_relu(W·[h v h u ])(9) in, The formula is (8) to calculate the feature vector of the heterogeneous edge, (9) to calculate the feature vector of the atom pair, (10) to calculate the weight of the neighbor atom to the target atom, (11) to calculate the context feature vector of the target atom, (12) to calculate the aggregation context feature vector of the target atom, (13) to calculate the latest state vector of the target atom, and (14) to predict the toxicity of the compound based on the latest state vector of the target atom.

Images