Semi-template drug molecule inverse synthesis prediction method based on graph-sequence pre-training

Abstract

The invention relates to a semi-template drug molecule inverse synthesis prediction method based on graph-sequence pre-training. The method comprises the following steps: selecting a public chemical molecule data set N; selecting a molecular sample from the N, and then extracting a molecular graph and a word segmentation sequence from the molecular sample; obtaining a sequence embedding matrix and a graph embedding matrix of the i by utilizing the word segmentation sequence of the i and the molecular graph; constructing an overall model M and a loss function Ltotal used for training the M; and training M in two stages by using AdamW to obtain a trained model M. By adopting the method disclosed by the invention, reactant atoms for generating a product can be accurately reduced, a starting material of the product is defined, and the research and development success rate of medicines or other products is improved.

Description

Technical Field

[0001] The present invention relates to the field of computer-aided drug synthesis, and in particular to a semi-template drug molecule reverse synthesis prediction method based on graph-sequence pre-training. Background Art

[0002] In modern drug discovery, retrosynthetic prediction of drug molecules is a critical step in drug development. This involves gradually "disassembling" target molecules into precursors or intermediates until commercially available or easily prepared starting materials are identified. Traditional retrosynthetic prediction relies heavily on the experience and trial-and-error of biochemists, leading to challenges such as low efficiency and high costs in drug development. With the rapid advancement of computer technology, computer-aided synthesis planning (CASP) systems have improved the efficiency of retrosynthetic prediction and reduced the cost of trial-and-error. Early CASP systems relied primarily on manual coding by experts to construct reaction rule libraries. However, this manual coding approach has significant limitations: the accuracy and completeness of the reaction rules obtained are limited by the experts' mastery of domain knowledge. Furthermore, as the number of chemical reactions increases exponentially, the cost and difficulty of manually collecting rules increase dramatically. These CASP systems suffer from limitations in accuracy, efficiency, and flexibility, making them difficult to adapt to the increasingly complex R&D needs of real-world applications.

[0003] In recent years, with the advancement of deep learning, deep learning models have been used to analyze vast amounts of chemical reaction data, autonomously mining universal molecular features and chemical reaction rules, and automatically predicting reactant structures and feasible synthetic routes. These models can be broadly categorized as template-based, template-free, and semi-template approaches. Template-based approaches rely on expert-constructed reaction template libraries, identifying reaction centers and bond-breaking patterns through template matching during the prediction process. While these approaches offer good interpretability, their performance is highly dependent on the coverage and quality of the templates, often underperforming for novel or unseen reaction types. Template-free approaches, on the other hand, rely entirely on no prior templates and typically employ graph neural networks or sequence generation models to infer reactant structures directly from product molecules. These approaches offer improved generalization capabilities, but suffer from poor interpretability and difficulty ensuring reaction plausibility. Semi-template approaches combine the advantages of both template-based and template-free approaches: leveraging the template concept, they identify the reaction center as auxiliary information to guide bond-breaking in the target molecule to generate synthons. Template-free approaches then complete these synthons to generate reactants, retaining the interpretability of template-based approaches while also possessing the generalization capabilities of template-free approaches. Consequently, they exhibit excellent overall performance and promising application prospects.

[0004] However, existing retrosynthesis prediction methods still face the following problems:

[0005] (1) Insufficient utilization of multidimensional information of drug molecules and limited molecular representation capabilities: Drug molecules can be represented by one-dimensional SMILES sequences and two-dimensional molecular graphs. Although SMILES sequence encoding captures linear grammatical structure, it omits molecular structural information such as chemical bond connections. Although molecular graph modeling can analyze structural information such as topological connectivity, it ignores sequence information such as atomic arrangement. Retrosynthesis prediction methods based on single-modal modeling cannot simultaneously take into account topological structure and sequence semantic information, resulting in the loss of chemical information and affecting the model's comprehensive learning of common molecular features.

[0006] (2) Scarcity of chemical reaction data and limited model generalization ability: Retrosynthesis prediction models are often trained using the finely calibrated open-source USPTO-50k dataset, but it only covers 50,000 reaction instances and is difficult to fully characterize more than 10 60 The chemical space of potential organic molecules constrains the model's ability to generalize to unseen molecules and rare reactions. Summary of the Invention

[0007] In response to the above-mentioned problems existing in the prior art, the technical problem to be solved by the present invention is: to solve the problem of limited molecular characterization and model generalization capabilities caused by insufficient utilization of multidimensional information of drug molecules and scarcity of reaction data, thereby improving the accuracy of retrosynthesis prediction of drug molecules.

[0008] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0009] A semi-template drug molecule retrosynthesis prediction method based on graph-sequence pre-training comprises the following steps:

[0010] S100: Select a public chemical molecule dataset, where the chemical molecule dataset contains N chemical molecules, each of which is a molecular sample and each molecular sample contains a SMILES string sequence;

[0011] S200: Select a molecular sample i from N, use a molecular parsing tool to extract atomic features and bond features from the SMILES string sequence in i, and then construct a molecular graph of i based on the atomic features and bond features; wherein the atomic features include the labeled real atom types, and the bond features include the labeled real bond types;

[0012] Perform word segmentation on the SMILES string sequence to obtain the word segmentation sequence of i;

[0013] Define the SMILES string sequence of i as the sequence mode of i, and define the molecular graph of i as the graph mode of i;

[0014] S300: Using the word segmentation sequence of i, obtain the sequence embedding matrix of i, and using the molecular graph of i, obtain the graph embedding matrix of i. The steps are as follows:

[0015] S310: Set the mask coverage ratios of the word segmentation sequence and the molecular graph of i respectively. According to the mask coverage ratios, the molecular graph of i is masked several times to obtain the masked graph enhanced view of i. The word segmentation sequence of i is masked several times to obtain the masked word segmentation sequence of i with ID. Embedding is used to process the masked word segmentation sequence of i with ID to obtain the initial sequence embedding representation of i. The masked areas covered by the masked graph enhanced view of i and the masked word segmentation sequence of i with ID do not overlap.

[0016] S320: Encode the image mask enhanced view of i and the initial sequence embedding representation of i using a graph encoder and a sequence encoder respectively: the image mask enhanced view of i is encoded by the graph encoder to obtain the image embedding matrix of i, and the initial sequence embedding representation of i is encoded by the sequence encoder to obtain the sequence embedding matrix of i;

[0017] Traverse all chemical molecules in N and repeat S200-S320 to obtain the graph embedding matrix and sequence embedding matrix corresponding to each molecule. The graph embedding matrix and sequence embedding matrix of all molecules constitute the pre-training set;

[0018] S330: Construct pre-training loss function L pretrain , take the pre-training set as input, use AdamW optimizer to pre-train the graph encoder and sequence encoder, use gradient descent method to reversely update the parameters of the graph encoder and sequence encoder, when L pretrain The training is stopped when the minimum value is reached or the training reaches the maximum number of iterations. At this time, the pre-trained graph encoder and pre-trained sequence encoder are obtained;

[0019] S400: Construct a retrosynthesis module, which includes a reaction center identification decoder and a synthon completion decoder; wherein the reaction center identification decoder includes an edge classifier, an atom classifier and a sigmoid function, which is used to construct a synthetic subgraph G sub The edge classifier and atom classifier are both double-layer MLPs; the synthon completion decoder is a Transformer-based graph editing decoder used to edit and generate molecular structure prediction graph G t ;

[0020] S500: The molecular graph of i obtained in S200 is used as the input of the retrosynthesis module and sequentially passes through the reaction center recognition decoder and synthon completion decoder to output the molecular structure prediction graph G of i. t ;

[0021] S600: Construct classification cross entropy loss L from edge vectors and atomic vectors center and edit sequence loss L edit, used to train the overall model M, which includes the pre-trained graph encoder and inverse synthesis module, as follows:

[0022] S610: Construct the total loss function L of M total , the calculation expression is as follows:

[0023] L total =λ1L center +λ2L edit

[0024] L center =λ bond ·BCE(p bond ,y bond )+λ atom ·BCE(p atom ,y atom )

[0025]

[0026] Among them, p bond represents the reaction probability prediction of each edge vector, p atom represents the reaction probability prediction for each atom vector, y bond and y atom Indicates whether the edge vector and the atomic vector are true labels, λ bond and λ atom represents an adjustable hyperparameter, BCE(·) represents the binary cross entropy loss function in the two-classification task, CE(·) represents the standard cross entropy loss function in the multi-classification task, represents the probability distribution of the edit operation target predicted at step t, represents the real node label or real key label corresponding to the target in step t, represents the probability distribution of the graph editing operation type predicted at step t, represents the actual graph editing operation type label at step t, represents the probability distribution of the operating parameters predicted at step t, represents the true parameter label corresponding to the t-th step, t represents the number of editing steps, T represents the maximum number of editing steps in synthesizer completion; λ1 and λ2 represent adjustable hyperparameters;

[0027] S720: Using L total And use AdamW to train M in two stages, as follows:

[0028] Phase 1: Set the total number of epochs, freeze the pre-trained graph encoder parameters, and select the first D epochs to train only the inverse synthesis module;

[0029] Phase 2: Starting from epoch D+1, a layer of the pre-trained graph encoder is unfrozen every fixed number of epochs, and then jointly trained with the inverse synthesis module that has been trained for D rounds. Gradient update permissions are gradually released in the pre-trained graph encoder layer order from bottom to top. When all pre-trained graph encoder layers are unfrozen and the total number of epochs is completed, training is stopped, and the trained overall model M' is obtained.

[0030] S800: Select the product molecule to be predicted and input M', and output the molecular structure prediction graph G of the product molecule to be predicted t ″,G t "A diagram of the molecules of reactants used to produce a product.

[0031] Preferably, in S200, the step of obtaining the molecular graph of i is as follows:

[0032] Use RDKit to parse the SMILES of i, convert the SMILES string into a molecular object using the Chem.MolFromSmiles function, and extract the atomic and bond features from the molecular object;

[0033] Map atomic features to node vectors, bond features to edge vectors, and the node vectors and edge vectors constitute the molecular graph of i.

[0034] Preferably, the content of the word segmentation sequence i generated by performing the word segmentation operation on the SMILES string sequence in S200 is as follows:

[0035] For the SMILES string sequence of i, the BPE word segmentation algorithm is used to construct several word lists, each of which contains several tokens. Each token is processed by the torch.nn.Embedding function to finally obtain the word segmentation sequence of i.

[0036] Preferably, in S320, the graph encoder is a GMPN graph encoder; and the sequence encoder is a Transformer encoder.

[0037] As an advantage, the pre-training loss function L is constructed in S330. pretrain The content is as follows:

[0038] S331: Use two layers of MLP to process node vectors and edge vectors, and construct a graph modality reconstruction loss function L based on the processing results graph-mask , L graph-mask The calculation expression is as follows:

[0039]

[0040] Among them, M gRepresents the set of masked node vectors in the molecular graph, v represents the node vector, Represents a set of adjacent nodes, u represents the adjacent node of v, a v represents the real atom corresponding to v, e u,v represents the real bond between v and u, h v represents the embedding representation of v, h u,v represents the embedding representation of the edge between u and v, and p(·) represents the conditional probability function;

[0041] Embed h for each mask position in the molecular graph of i using a linear projection layer t Perform classification and construct the sequence modality reconstruction loss function L according to the classification results seq-mask , L seq-mask The calculation expression is as follows:

[0042]

[0043] in, Represents the masked word vector set in the word segmentation sequence, c represents the masked word vector, w c represents the real atom corresponding to c, h c represents the embedded representation of c;

[0044] S332: Construct the total loss function L of the graph mode and sequence mode of i recon , the calculation expression is as follows:

[0045] L recon =λ g L graph-mask +λ s L seq-mask

[0046] Among them, λ g ,λ s Both represent hyperparameters, which are used to control the weights of their corresponding modalities;

[0047] S333: Using InfoNCE to construct a cross-modal graph loss function L gg and the cross-modal sequence loss function L ss :

[0048] Set the mapping vector dimension value d z , the graph modality of i and the sequence modality of i are respectively input into their corresponding two-layer MLP for nonlinear mapping, and both are activated by LeakyReLU activation function, and then the corresponding embedding vectors are output: the graph modality of i corresponds to the output graph embedding vector The sequence mode of i corresponds to the output sequence embedding vector L gg The calculation expression is as follows:

[0049]

[0050] in, represents the set graph query vector, represents the key vector of the graph positive sample, represents the graph negative sample bond vector, τ represents the temperature coefficient; the positive sample is: select two graph views in molecule i and when When is the query vector, The positive sample is defined as i; the negative sample is: randomly select a molecule j other than i in N, any graph view in j It is considered to be a negative sample of i;

[0051] Similarly, L ss The calculation expression is as follows:

[0052]

[0053] in, represents the sequence query vector, Represents the sequence positive sample key vector, z sj′ (j′≠i′,j′∈N) represents the graph negative sample key vector;

[0054] InfoNCE loss is used to construct a cross-modal graph sequence loss function L corss The calculation expression is as follows:

[0055]

[0056] Calculate the total cross-modal comparison loss function L contrast , the calculation expression is as follows:

[0057] L contrast =λ gg L gg +λ ss L ss +λ cross L cross

[0058] Among them, λ gg ,λ ss ,λ cross Both represent the loss weights corresponding to each loss function;

[0059] S334: Joint L recon and L contrast Get the pre-training loss function L pretrain , the calculation formula is as follows:

[0060] Lpretrain =L recon +L contrast .

[0061] As an example, the molecular structure prediction graph G of i is output in S500. t The steps are as follows:

[0062] S510: Input all edge vectors in the molecular graph of i into the edge classifier, and output the probability p of each edge vector being a broken bond bond (e) Input all atomic features in the molecular graph of i into the atom classifier and output the probability p of the atom corresponding to the atomic feature participating in the reaction. atom (a); then the comprehensive breaking score of the edge between the two atoms is Score(e);

[0063] Select a u and a v , a u Indicates that the node u corresponds to a real atom, then a u and a v The edge between the two is a broken bond, and the calculation formula for the comprehensive fracture score is as follows:

[0064] Score(e)=p bond (e)·p atom (au)·p atom (a v )

[0065] Where, e=(a u ,a v ) means a u with a v The edge between

[0066] S520: Repeat S510 to calculate the comprehensive fracture scores of all edges in i, sort the comprehensive fracture scores of all edges in descending order, and select the edges corresponding to the first n comprehensive fracture scores as the fracture keys of i;

[0067] S530: break n broken bonds, and the associated atoms at both ends of each broken bond and the first-order neighboring structures corresponding to the two associated atoms constitute the synthetic subgraph G after each broken bond is broken. sub ;

[0068] S540: Each G sub Input the synthesizer completion decoder to complete the synthesizer graph and obtain the current molecular structure prediction graph G of i t .

[0069] As a preference, in the step S540, each G sub The content input to the synthesizer completion decoder for synthesizer subgraph completion is as follows:

[0070] S541: Define the graph editing operation dictionary structure A. The expression of A is as follows:

[0071] A={AddNode,RemoveNode,AddBond,ChangeBond,AddGroup,Terminate}

[0072] Among them, AddNond indicates atom addition, removeNode indicates atom deletion, AddBond indicates bond formation, ChangeBond indicates bond breaking, AddGroup indicates group addition, and Terminate indicates termination;

[0073] S542: Each G sub Input synthesizer completion decoder, using A to G sub Perform graph editing and obtain the editing result EditOp, which is expressed as follows:

[0074] EditOp=<target,op,param>

[0075] Where target represents the node or edge to be edited, op represents the operation type of graph editing, and param represents the newly added atom type, bond type, or group type.

[0076] S543: Update the input G according to EditOp sub Get G′ sub , and for G′ sub Perform conditional judgment, the judgment conditions are as follows:

[0077] 1) All currently broken atoms are saturated and there are no isolated nodes; 2) The molecule has no charge anomalies or valence conflicts; 3) The cumulative number of edit operations on a single synthetic subgraph exceeds T;

[0078] When any of the above three conditions is met, the graph editing operation is stopped, the completion operation of the composite subgraph is completed, and G′ sub As molecular structure prediction graph G t Output, otherwise G′ sub =G sub And return to S542.

[0079] Compared with the prior art, the present invention has at least the following advantages:

[0080] 1. This paper simultaneously utilizes two modalities, molecular graph structure and SMILES sequence, and collaboratively models them through graph neural networks and Transformer encoders. It also introduces intra-modal reconstruction tasks to guide the model to fully understand the semantics of chemical molecules from both the structural and linguistic levels, achieving information complementarity and structure-semantic synergy between graphs and sequences, overcoming the limited expressive power of existing single-modal methods.

[0081] 2. A multi-view contrastive learning mechanism is proposed. Random masks are introduced in training to generate multiple views. Three types of positive sample pairs are constructed: image-image, sequence-sequence, and image-sequence. Contrastive learning is used to maximize the representation consistency of the same molecule in different views for modal alignment, effectively alleviating the problem of training instability in scenarios with small samples or low-quality reaction data.

[0082] 3. During the retrosynthetic generation process, the present invention adopts a semi-templated retrosynthesis method, explicitly introduces a reaction center identification module, integrates edge and atom classification information to accurately predict the break position, and constructs an initial synthetic subgraph based on the break center. Furthermore, by defining a graph editing operation space, autoregressive generation and beam search strategies are used to gradually complete the synthetic substructure to form reactants. This process improves the interpretability and chemical rationality of the reaction path.

[0083] 4. Through a freeze-thaw parameter scheduling strategy, the pre-trained encoder is frozen in the first stage, and only the decoder and reaction center identification module are trained to avoid premature interference with the learned graph-sequence representation. In the second stage, the encoder is unfrozen layer by layer and jointly optimized with the decoder to promote a more comprehensive transfer of pre-trained knowledge to the retrosynthesis task, achieving effective coordination between the two sub-tasks of reaction center identification and reactant editing. This can not only accelerate convergence but also prevent gradient interference between tasks, with significant advantages in training efficiency and multi-task fusion effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 Schematic diagram of the method of the present invention.

[0085] Figure 2 Schematic diagram of the molecular characterization method of the present invention. DETAILED DESCRIPTION

[0086] The present invention is described in further detail below.

[0087] See also Figure 1-Figure 2 A semi-template drug molecule retrosynthesis prediction method based on graph-sequence pre-training includes the following steps:

[0088] S100: Select a public chemical molecule dataset, where the chemical molecule dataset contains N chemical molecules, each of which is a molecular sample and each molecular sample contains a SMILES string sequence;

[0089] S200: Select a molecular sample i from N, use a molecular parsing tool to extract atomic features and bond features from the SMILES string sequence in i, and then construct a molecular graph of i based on the atomic features and bond features; wherein the atomic features include the labeled real atom types, and the bond features include the labeled real bond types;

[0090] Perform a word segmentation operation on the SMILES string sequence to obtain a word segmentation sequence of i; the word segmentation operation is a prior art;

[0091] Define the SMILES string sequence of i as the sequence mode of i, and define the molecular graph of i as the graph mode of i;

[0092] The content of the word segmentation sequence i generated by performing word segmentation operation on the SMILES string sequence in S200 is as follows:

[0093] For the SMILES string sequence of i, the BPE word segmentation algorithm is used to construct several word lists, each of which contains several tokens. Each token is processed by the torch.nn.Embedding function to finally obtain the word segmentation sequence of i.

[0094] To unify the model input, all SMILES sequences are padded or truncated to a maximum length of 200. At the same time, SMILES string sequences that are chemically invalid or fail to parse are filtered to improve data quality. The BPE word segmentation algorithm and torch.nn.Embedding are both existing technologies.

[0095] In S200, the steps of obtaining the molecular graph of i are as follows:

[0096] Use RDKit to parse the SMILES of i, convert the SMILES string into a molecular object using the Chem.MolFromSmiles function, and extract the atomic and bond features from the molecular object;

[0097] Map atomic features to node vectors, bond features to edge vectors, and the node vectors and edge vectors constitute the molecular graph of i.

[0098] Both RDKit and Chem.MolFromSmiles functions are existing technologies; atomic features include atom type, formal charge, hybridization state, valence state, number of hydrogen atoms, and chirality information; bond features include bond type, whether it is conjugated, whether it is aromatic, and whether it is cyclic.

[0099] S300: Using the word segmentation sequence of i, obtain the sequence embedding matrix of i, and using the molecular graph of i, obtain the graph embedding matrix of i. The steps are as follows:

[0100] S310: Set the mask coverage ratios of the word segmentation sequence and the molecular graph of i respectively. According to the mask coverage ratios, the molecular graph of i is masked several times to obtain the masked graph enhanced view of i. The word segmentation sequence of i is masked several times to obtain the masked word segmentation sequence of i with ID. Embedding is used to process the masked word segmentation sequence of i with ID to obtain the initial sequence embedding representation of i. The masked areas covered by the masked graph enhanced view of i and the masked word segmentation sequence of i with ID do not overlap.

[0101] Generally, for graph modality, we use NetworkX to calculate node degrees and generate weighted sampling probabilities, prioritizing masking high-connectivity nodes and setting their node and edge features to zero. For sequence modality, we count the token frequencies in the word segmentation sequence, preferentially sample 15% of low-frequency tokens and replace them with [MASK] to ensure that the masked areas of the two modalities do not overlap.

[0102] S320: Encode the image mask enhanced view of i and the initial sequence embedding representation of i using a graph encoder and a sequence encoder respectively: the image mask enhanced view of i is encoded by the graph encoder to obtain the image embedding matrix of i, and the initial sequence embedding representation of i is encoded by the sequence encoder to obtain the sequence embedding matrix of i;

[0103] Traverse all chemical molecules in N and repeat S200-S320 to obtain the graph embedding matrix and sequence embedding matrix corresponding to each molecule. The graph embedding matrix and sequence embedding matrix of all molecules constitute the pre-training set;

[0104] In S320, the graph encoder is a GMPN graph encoder; the sequence encoder is a Transformer encoder; both the GMPN graph encoder and the Transformer encoder are prior arts.

[0105] S330: Construct pre-training loss function L pretrain , take the pre-training set as input, use AdamW optimizer to pre-train the graph encoder and sequence encoder, use gradient descent method to reversely update the parameters of the graph encoder and sequence encoder, when L pretinTraining is stopped when the minimum value is reached or the maximum number of iterations has been reached, at which point a pre-trained graph encoder and a pre-trained sequence encoder are obtained. The AdamW optimizer is a state-of-the-art technique. During the learning process, the graph encoder is forced to adjust its internal parameters and representation method so that the graph embedding matrix information it generates can effectively "predict" or "match" the sequence embedding matrix information of the same molecule generated by the sequence encoder. This process is equivalent to using sequence information as a supervisory signal, guiding the graph encoder to encode in its generated graph embedding those features that are useful for distinguishing molecular identities, and these features are precisely those captured by the sequence encoder (such as functional group type, approximate connection order, molecular size, certain key chemical substructures, etc.).

[0106] In S330, a pre-training loss function L is constructed. pretrain The content is as follows:

[0107] S331: Use two layers of MLP to process node vectors and edge vectors, and construct a graph modality reconstruction loss function L based on the processing results graph-mask , L graph-mask The calculation expression is as follows:

[0108]

[0109] Among them, M g Represents the set of masked node vectors in the molecular graph, v represents the node vector, Represents a set of adjacent nodes, u represents the adjacent node of v, a v represents the real atom corresponding to v, e u,v represents the real bond between v and u, h v represents the embedding representation of v, h u,v represents the embedding representation of the edge between u and v, p(·) represents the conditional probability function; the two-layer MLP represents the two-layer multi-head attention mechanism, which is the existing technology;

[0110] Embed h for each mask position in the molecular graph of i using a linear projection layer t Perform classification and construct the sequence modality reconstruction loss function L according to the classification results seq-mask , L seq-mask The calculation expression is as follows:

[0111]

[0112] in, Represents the masked word vector set in the word segmentation sequence, c represents the masked word vector, w c represents the real atom corresponding to c, h c represents the embedding representation of c; the linear projection layer is the existing technology;

[0113] S332: Construct the total loss function L of the graph mode and sequence mode of i recon , the calculation expression is as follows:

[0114] L recon =λ g L graph-mask +λ s L seq-mask

[0115] Among them, λ g ,λ s Both represent hyperparameters, which are used to control the weights of their corresponding modalities;

[0116] S333: Using InfoNCE to construct a cross-modal graph loss function L gg and the cross-modal sequence loss function L ss , contrast loss InfoNCE is the existing technology:

[0117] Set the mapping vector dimension value d z , the graph modality of i and the sequence modality of i are respectively input into their corresponding two-layer MLP for nonlinear mapping, and both are activated by LeakyReLU activation function, and then the corresponding embedding vectors are output: the graph modality of i corresponds to the output graph embedding vector The sequence mode of i corresponds to the output sequence embedding vector L gg The calculation expression is as follows:

[0118]

[0119] in, represents the set graph query vector, represents the key vector of the graph positive sample, represents the graph negative sample bond vector, τ represents the temperature coefficient; the positive sample is: select two graph views in molecule i and when When is the query vector, The positive sample is defined as i; the negative sample is: randomly select a molecule j other than i in N, any graph view in j It is considered to be a negative sample of i;

[0120] Similarly, L ss The calculation expression is as follows:

[0121]

[0122] in, represents the sequence query vector, Represents the sequence positive sample key vector, z sj′(j′≠i′,j′∈N) represents the graph negative sample key vector;

[0123] InfoNCE loss is used to construct a cross-modal graph sequence loss function L cross The calculation expression is as follows:

[0124]

[0125] Calculate the total cross-modal comparison loss function L contrast , the calculation expression is as follows:

[0126] L contrast =λ gg L gg +λ ss L ss +λ cross L cross

[0127] Among them, λ gg ,λ ss ,λ cross Both represent the loss weights corresponding to each loss function;

[0128] S334: Joint L recon and L contras tGet the pre-training loss function L pretrain , the calculation formula is as follows:

[0129] L pretrain =L recon +L contrast .

[0130] The core principle of pre-training is to use the alignment effect of contrastive learning: the goal in contrastive learning is to make the graph embedding and sequence embedding of the same molecule (positive sample) as close as possible in the representation space (maximize similarity), while the embeddings of different molecules (negative samples) are as far apart as possible (minimize similarity); to achieve this goal, the graph encoder is forced to adjust its internal parameters and representation methods during the learning process so that the graph embedding it produces can effectively "predict" or "match" the sequence embedding of the same molecule produced by the sequence encoder; this process is equivalent to using sequence information as a supervisory signal, guiding the graph encoder to encode those features that are useful for distinguishing the identity of the molecule in its generated graph embedding, which are precisely the features that the sequence encoder also captures (such as functional group type, approximate connection order, molecular size, certain key chemical substructures, etc.); successful contrastive learning forces the graph embedding and sequence embedding to express the molecule in the same shared semantic space. The output of the graph encoder is no longer just an abstraction of the graph structure, but incorporates an understanding of the chemical nature of the molecule, which is partly obtained through alignment with the sequence information.

[0131] S400: Construct a retrosynthesis module, which includes a reaction center identification decoder and a synthon completion decoder; wherein the reaction center identification decoder includes an edge classifier, an atom classifier and a sigmoid function, which is used to construct a synthetic subgraph G sub The edge classifier and atom classifier are both double-layer MLPs; the synthon completion decoder is a Transformer-based graph editing decoder used to edit and generate molecular structure prediction graph G t ; MLP is a multi-layer perceptron and Transformer are both existing technologies;

[0132] S500: The molecular graph of i obtained in S200 is used as the input of the retrosynthesis module and sequentially passes through the reaction center recognition decoder and synthon completion decoder to output the molecular structure prediction graph G of i. t ;

[0133] The molecular structure prediction graph G of i is output in S500 t The steps are as follows:

[0134] S510: Input all edge vectors in the molecular graph of i into the edge classifier, and output the probability p of each edge vector being a broken bond bond (e) Input all atomic features in the molecular graph of i into the atom classifier and output the probability p of the atom corresponding to the atomic feature participating in the reaction. atom (a); then the comprehensive breaking score of the edge between the two atoms is Score(e);

[0135] Select a u and a v , a u Indicates that the node u corresponds to a real atom, then a u and a v The edge between the two is a broken bond, and the calculation formula for the comprehensive fracture score is as follows:

[0136] Score(e)=p bond (e)·p atom (a u )·p atom (a v )

[0137] Where, e=(a u ,a v ) means a u with a v The edge between

[0138] S520: Repeat S510 to calculate the comprehensive fracture scores of all edges in i, sort the comprehensive fracture scores of all edges in descending order, and select the edges corresponding to the first n comprehensive fracture scores as the fracture keys of i;

[0139] S530: break n broken bonds, and the associated atoms at both ends of each broken bond and the first-order neighboring structures corresponding to the two associated atoms constitute the synthetic subgraph G after each broken bond is broken. sub ;

[0140] S540: Each G sub Input the synthesizer completion decoder to complete the synthesizer graph and obtain the current molecular structure prediction graph G of i t ;

[0141] The S540 will each G sub The content input to the synthesizer completion decoder for synthesizer subgraph completion is as follows:

[0142] S541: Define the graph editing operation dictionary structure A. The expression of A is as follows:

[0143] A={AddNode,RemoveNode,AddBond,ChangeBond,AddGroup,Terminate}

[0144] Among them, AddNode indicates atom addition, removeNode indicates atom deletion, AddBond indicates bond formation, ChangeBond indicates bond breaking, AddGroup indicates group addition, and Terminate indicates termination;

[0145] S542: Each G sub Input synthesizer completion decoder, using A to G sub Perform graph editing and obtain the editing result EditOp, which is expressed as follows:

[0146] EditOp=<target,op,param>

[0147] Where target represents the node or edge to be edited, op represents the operation type of graph editing, and param represents the newly added atom type, bond type, or group type.

[0148] S543: Update the input G according to EditOp sub Get G′ sub , and for G′ sub Perform conditional judgment, the judgment conditions are as follows:

[0149] 1) All currently broken atoms are saturated and there are no isolated nodes; 2) The molecule has no charge anomalies or valence conflicts; 3) The cumulative number of edit operations on a single synthetic subgraph exceeds T;

[0150] When any of the above three conditions is met, the graph editing operation is stopped, the completion operation of the composite subgraph is completed, and G′ sub As molecular structure prediction graph G t Output, otherwise G′ sub =G sub And return to S542.

[0151] S600: Construct classification cross entropy loss L from edge vectors and atomic vectors center and edit sequence loss L edit , used to train the overall model M, which includes a pre-trained graph encoder and an inverse synthesis module. Since the graph encoder has been optimized in the pre-training stage, its graph embedding can well "represent" the corresponding sequence embedding. Therefore, in the fine-tuning stage, even if the sequence encoder is not used, the graph embedding generated by the graph encoder itself implies (part of) the information learned by the sequence encoder. This information is "distilled" into the parameters of the graph encoder and its output embedding vector, as follows:

[0152] S610: Construct the total loss function L of M total , the calculation expression is as follows:

[0153] L total =λ1L center +λ2L edit

[0154] L center =λ bond ·BCE(p bond ,y bond )+λ atom ·BCE(p atom ,y atom )

[0155]

[0156] Among them, p bond represents the reaction probability prediction of each edge vector, p atom represents the reaction probability prediction for each atom vector, y bond and y atom Indicates whether the edge vector and the atomic vector are real labels, 0 means false, 1 means true, λ bond and λ atom represents an adjustable hyperparameter, BCE(·) represents the binary cross entropy loss function in the two-classification task, CE(·) represents the standard cross entropy loss function in the multi-classification task, represents the probability distribution of the edit operation target predicted at step t, represents the real node label or real key label corresponding to the target in step t, represents the probability distribution of the graph editing operation type predicted at step t, represents the actual graph editing operation type label at step t, represents the probability distribution of the operating parameters predicted at step t, represents the true parameter label corresponding to the t-th step, t represents the number of editing steps, T represents the maximum number of editing steps in synthesizer completion; λ1 and λ2 represent adjustable hyperparameters;

[0157] BCE(·) represents the binary cross entropy loss function in the binary classification task, CE(·) represents the standard cross entropy loss function in the multi-classification task, represents the target (atom or bond) probability distribution of the edit operation predicted at step t, represents the real target node or key label corresponding to the t-th step, represents the probability distribution of the graph editing operation type (atom addition, atom deletion, bond formation, bond breaking, group addition, termination) predicted at step t, Indicates the actual operation type label corresponding to step t, represents the probability distribution of the operational parameters (new atom, bond, or group type) predicted or updated at step t, represents the true parameter label corresponding to the t-th step, t represents the number of editing steps, t∈T, T represents the maximum number of editing steps in the synthesizer completion;

[0158] S720: Using L total And use AdamW to train M in two stages. AdamW is an existing technology, and the content is as follows:

[0159] Phase 1: Set the total number of epochs, freeze the pre-trained graph encoder parameters, and select the first D epochs to train only the retrosynthesis module. The purpose of this phase is to prevent catastrophic forgetting: the pre-trained graph encoder has already learned common molecular structures and semantic representations in large-scale unsupervised or multimodal tasks. If small-scale supervised tasks (such as reaction center identification) are used directly in the early stages of fine-tuning, the full parameter linkage training will result in training failing to achieve the desired results. Therefore, the encoder parameters should be frozen initially, and only the decoder in the retrosynthesis module should be trained to allow it to first learn how to perform tasks based on the existing representations without "contaminating" the pre-trained training capabilities.

[0160] Phase 2: Starting from epoch D+1, a layer of the pre-trained graph encoder is unfrozen every fixed number of rounds, and the inverse synthesis module that has been trained for D rounds is jointly trained. Gradient update rights are gradually released in the order of layers from bottom to top in the pre-trained graph encoder. When all pre-trained graph encoder layers are unfrozen and the total number of epochs is completed, training is stopped, and the trained overall model M' is obtained. The purpose of unfreezing layer by layer is to allow the encoder to gradually adapt to downstream tasks while retaining the pre-training ability. Instead of freezing it forever, the parameter update rights are gradually released layer by layer, so that the model can stably transfer existing high-quality information. Avoiding conflicts between new and old parameters during gradient updates and gradually adapting to fine-tuning results to enhance model adaptability is a typical technique of "progressive transfer learning," which has been widely proven effective in natural language processing (NLP), computer vision, and graph learning. This part of the training can be considered the fine-tuning phase of the model: only the graph encoder from pre-training is used, because its results contain implicit sequence information. Since the graph encoder has been optimized during pre-training, its graph embeddings can well "represent" the corresponding sequence embeddings. Therefore, during the fine-tuning phase, even without using the sequence encoder, the graph embeddings generated by the graph encoder implicitly contain (partial) information learned by the sequence encoder. This information is "distilled" into the graph encoder parameters and its output embedding vector.

[0161] S800: Select the product molecule to be predicted and input M', and output the molecular structure prediction graph G of the product molecule to be predicted t ″,G t "A diagram of the molecules of reactants used to produce a product.

[0162] Experimental content and results

[0163] Experiments were conducted on the standard public dataset USPTO-50k, which contains approximately 50,000 labeled reaction data and is divided into training, validation, and test sets in a ratio of 40,000 / 5,000 / 5,000. The SMILES sequence was parsed into a molecular graph using the RDKit tool, and multidimensional graph structure features such as atom type, hybridization state, and bond type were extracted. The pre-trained GMPN encoder was then migrated to a semi-template model architecture, using a freeze-thaw joint training strategy: the first four encoder layers were frozen for the first 10 training cycles, and only the decoder was optimized; thereafter, one layer was thawed every five cycles, and parameter updates were gradually released.

[0164] The top-k accuracy (k = 1, 3, 5, 10) is used as the primary evaluation metric. This is defined as whether the top k candidate reactants output by the model contain the actual reactant structure. This metric is used to measure the reliability and usability of the model's predictions.

[0165] In summary, the present invention proposes a semi-template drug molecule retrosynthesis prediction method based on graph-sequence pre-training, implements non-overlapping masking and contrastive learning to fuse molecular graphs and sequences for molecular representation pre-training, and performs two-stage fine-tuning optimization and performance evaluation based on reaction center identification and synthon completion tasks. On the one hand, it solves the problem of limited molecular representation caused by insufficient utilization of multi-dimensional information of drug molecules, and on the other hand, it solves the problem of limited model generalization ability caused by scarce reaction data, and finally realizes drug molecule retrosynthesis prediction with high accuracy, good generalization ability and interpretability. The present invention can be applied to actual drug molecule retrosynthesis prediction to accelerate the drug research and development process.

[0166] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the purpose and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A semi-template drug molecule retrosynthesis prediction method based on graph-sequence pre-training, characterized by: The steps include: S100: Select a public chemical molecule dataset, where the chemical molecule dataset contains N chemical molecules, each of which is a molecular sample and each molecular sample contains a SMILES string sequence; S200: Select a molecular sample i from N, use a molecular parsing tool to extract atomic features and bond features from the SMILES string sequence in i, and then construct a molecular graph of i based on the atomic features and bond features; wherein the atomic features include the labeled real atom types, and the bond features include the labeled real bond types; Perform word segmentation on the SMILES string sequence to obtain the word segmentation sequence of i; Define the SMILES string sequence of i as the sequence mode of i, and define the molecular graph of i as the graph mode of i; S300: Using the word segmentation sequence of i, obtain the sequence embedding matrix of i, and using the molecular graph of i, obtain the graph embedding matrix of i. The steps are as follows: S310: Set the mask coverage ratios of the word segmentation sequence and the molecular graph of i respectively. According to the mask coverage ratios, the molecular graph of i is masked several times to obtain the masked graph enhanced view of i. The word segmentation sequence of i is masked several times to obtain the masked word segmentation sequence of i with ID. Embedding is used to process the masked word segmentation sequence of i with ID to obtain the initial sequence embedding representation of i. The masked areas covered by the masked graph enhanced view of i and the masked word segmentation sequence of i with ID do not overlap. S320: Encode the image mask enhanced view of i and the initial sequence embedding representation of i using a graph encoder and a sequence encoder respectively: the image mask enhanced view of i is encoded by the graph encoder to obtain the image embedding matrix of i, and the initial sequence embedding representation of i is encoded by the sequence encoder to obtain the sequence embedding matrix of i; Traverse all chemical molecules in N and repeat S200-S320 to obtain the graph embedding matrix and sequence embedding matrix corresponding to each molecule. The graph embedding matrix and sequence embedding matrix of all molecules constitute the pre-training set; S330: Construct pre-training loss function L pretrain , take the pre-training set as input, use AdamW optimizer to pre-train the graph encoder and sequence encoder, use gradient descent method to reversely update the parameters of the graph encoder and sequence encoder, when L pretrain The training is stopped when the minimum value is reached or the training reaches the maximum number of iterations. At this time, the pre-trained graph encoder and pre-trained sequence encoder are obtained; S400: Construct a retrosynthesis module, which includes a reaction center identification decoder and a synthon completion decoder; wherein the reaction center identification decoder includes an edge classifier, an atom classifier and a sigmoid function, which is used to construct a synthetic subgraph G sub The edge classifier and atom classifier are both double-layer MLPs; the synthon completion decoder is a Transformer-based graph editing decoder used to edit and generate molecular structure prediction graph G t ; S500: The molecular graph of i obtained in S200 is used as the input of the retrosynthesis module and sequentially passes through the reaction center recognition decoder and synthon completion decoder to output the molecular structure prediction graph G of i. t ; S600: Construct classification cross entropy loss L from edge vectors and atomic vectors center and edit sequence loss L edit , used to train the overall model M, which includes the pre-trained graph encoder and inverse synthesis module, as follows: S610: Construct the total loss function L of M total , the calculation expression is as follows: L total =λ1L center +λ2L edit L center =λ bond· BCE(p bond ,y bond )+λ atom ·BCE(p atom ,y atom ) Among them, p bond represents the reaction probability prediction of each edge vector, p atom represents the reaction probability prediction for each atom vector, y bond and y atom Indicates whether the edge vector and the atomic vector are true labels, λ bond and λ atom represents an adjustable hyperparameter, BCE(·) represents the binary cross entropy loss function in the two-classification task, CE(·) represents the standard cross entropy loss function in the multi-classification task, represents the probability distribution of the edit operation target predicted at step t, represents the real node label or real key label corresponding to the target in step t, represents the probability distribution of the graph editing operation type predicted at step t, represents the actual graph editing operation type label at step t, represents the probability distribution of the operating parameters predicted at step t, represents the true parameter label corresponding to the t-th step, t represents the number of editing steps, T represents the maximum number of editing steps in synthesizer completion; λ1 and λ2 represent adjustable hyperparameters; S720: Using L total And use AdamW to train M in two stages, as follows: Phase 1: Set the total number of epochs, freeze the pre-trained graph encoder parameters, and select the first D epochs to train only the inverse synthesis module; Phase 2: Starting from epoch D+1, a layer of the pre-trained graph encoder is unfrozen every fixed number of epochs, and then jointly trained with the inverse synthesis module that has been trained for D rounds. Gradient update permissions are gradually released in the pre-trained graph encoder layer order from bottom to top. When all pre-trained graph encoder layers are unfrozen and the total number of epochs is completed, training is stopped, and the trained overall model M' is obtained. S800: Select the product molecule to be predicted and input M', and output the molecular structure prediction graph G of the product molecule to be predicted t ″,G t "A diagram of the molecules of reactants used to produce a product.

2. The method for predicting the retrosynthesis of drug molecules using a semi-template method based on graph-sequence pre-training according to claim 1, wherein: In S200, the steps of obtaining the molecular graph of i are as follows: Use RDKit to parse the SMILES of i, convert the SMILES string into a molecular object using the Chem.MolFromSmiles function, and extract the atomic and bond features from the molecular object; Map atomic features to node vectors, bond features to edge vectors, and the node vectors and edge vectors constitute the molecular graph of i.

3. The method for predicting the retrosynthesis of a semi-templated drug molecule based on graph-sequence pre-training according to claim 2, wherein: The content of the word segmentation sequence i generated by performing word segmentation operation on the SMILES string sequence in S200 is as follows: For the SMILES string sequence of i, the BPE word segmentation algorithm is used to construct several word lists, each of which contains several tokens. Each token is processed by the torch.nn.Embedding function to finally obtain the word segmentation sequence of i.

4. The method for predicting the retrosynthesis of drug molecules using a semi-template method based on graph-sequence pre-training according to claim 3, wherein: In S320, the graph encoder is a GMPN graph encoder; and the sequence encoder is a Transformer encoder.

5. The method for predicting the retrosynthesis of drug molecules based on graph-sequence pre-training according to claim 4, characterized in that: In S330, a pre-training loss function L is constructed. pretrain The content is as follows: S331: Use two layers of MLP to process node vectors and edge vectors, and construct a graph modality reconstruction loss function L based on the processing results graph-mask , L graph-mask The calculation expression is as follows: Among them, M g Represents the set of masked node vectors in the molecular graph, v represents the node vector, Represents a set of adjacent nodes, u represents the adjacent node of v, a v represents the real atom corresponding to v, e u,v represents the real bond between v and u, h v represents the embedding representation of v, h u,v represents the embedding representation of the edge between u and v, and p(·) represents the conditional probability function; Embed h for each mask position in the molecular graph of i using a linear projection layer t Perform classification and construct the sequence modality reconstruction loss function L according to the classification results seq-mask , L seq-mask The calculation expression is as follows: in, Represents the masked word vector set in the word segmentation sequence, c represents the masked word vector, w c represents the real atom corresponding to c, h c represents the embedded representation of c; S332: Construct the total loss function L of the graph mode and sequence mode of i recon , the calculation expression is as follows: L recon =λ g L graph-mask +λ s L seq-mask Among them, λ g ,λ s Both represent hyperparameters, which are used to control the weights of their corresponding modalities; S333: Using InfoNCE to construct a cross-modal graph loss function L gg and the cross-modal sequence loss function L ss : Set the mapping vector dimension value d z , the graph modality of i and the sequence modality of i are respectively input into their corresponding two-layer MLP for nonlinear mapping, and both are activated by LeakyReLU activation function, and then the corresponding embedding vectors are output: the graph modality of i corresponds to the output graph embedding vector The sequence mode of i corresponds to the output sequence embedding vector L gg The calculation expression is as follows: in, represents the set graph query vector, represents the key vector of the graph positive sample, represents the graph negative sample bond vector, τ represents the temperature coefficient; the positive sample is: select two graph views in molecule i and when When is the query vector, The positive sample is defined as i; the negative sample is: randomly select a molecule j other than i in N, any graph view in j It is considered to be a negative sample of i; Similarly, L ss The calculation expression is as follows: in, represents the sequence query vector, Represents the sequence positive sample key vector, z sj′ (j′≠i′,j′∈N) represents the graph negative sample key vector; InfoNCE loss is used to construct a cross-modal graph sequence loss function L cross The calculation expression is as follows: Calculate the total cross-modal comparison loss function L contrast , the calculation expression is as follows: L contrast =λ gg L gg +λ ss L ss +λ cross L cross Among them, λ gg ,λ ss ,λ cross Both represent the loss weights corresponding to each loss function; S334: Joint L recon and L contrast Get the pre-training loss function L pretrain , the calculation formula is as follows: L pretrain =L recon +L contrast 。 6. The method for predicting the retrosynthesis of a semi-templated drug molecule based on graph-sequence pre-training according to claim 5, characterized in that: The molecular structure prediction graph G of i is output in S500 t The steps are as follows: S510: Input all edge vectors in the molecular graph of i into the edge classifier, and output the probability p of each edge vector being a broken bond bond (e) Input all atomic features in the molecular graph of i into the atom classifier and output the probability p of the atom corresponding to the atomic feature participating in the reaction. atom (a); then the comprehensive breaking score of the edge between the two atoms is Score(e); Select a u and a v , a u Indicates that the node u corresponds to a real atom, then a u and a v The edge between the two is a broken bond, and the calculation formula for the comprehensive fracture score is as follows: Score(e)=p bond (e)·p atom (a u )·p atom (a v ) Where, e=(a u ,a v ) means a u with a v The edge between S520: Repeat S510 to calculate the comprehensive fracture scores of all edges in i, sort the comprehensive fracture scores of all edges in descending order, and select the edges corresponding to the first n comprehensive fracture scores as the fracture keys of i; S530: break n broken bonds, and the associated atoms at both ends of each broken bond and the first-order neighboring structures corresponding to the two associated atoms constitute the synthetic subgraph G after each broken bond is broken. sub ; S540: Each G sub Input the synthesizer completion decoder to complete the synthesizer graph and obtain the current molecular structure prediction graph G of i t .

7. The method for predicting the retrosynthesis of a semi-templated drug molecule based on graph-sequence pre-training according to claim 6, characterized in that: The S540 will each G sub The content input to the synthesizer completion decoder for synthesizer subgraph completion is as follows: S541: Define the graph editing operation dictionary structure A. The expression of A is as follows: A={AddNode,RemoveNode,AddBond,ChangeBond,AddGroup,Terminate} Among them, AddNode indicates atom addition, removeNode indicates atom deletion, AddBond indicates bond formation, ChangeBond indicates bond breaking, AddGroup indicates group addition, and Terminate indicates termination; S542: Each G sub Input synthesizer completion decoder, using A to G sub Perform graph editing and obtain the editing result EditOp, which is expressed as follows: EditOp=<target,op,param> Where target represents the node or edge to be edited, op represents the operation type of graph editing, and param represents the newly added atom type, bond type, or group type. S543: Update the input G according to EditOp sub Get G′ sub , and for G′ sub Perform conditional judgment, the judgment conditions are as follows: 1) All currently broken atoms are saturated and there are no isolated nodes; 2) The molecule has no charge anomalies or valence conflicts; 3) The cumulative number of edit operations on a single synthetic subgraph exceeds T; When any of the above three conditions is met, the graph editing operation is stopped, the completion operation of the composite subgraph is completed, and G′ sub As molecular structure prediction graph G t Output, otherwise G′ sub =G sub And return to S542.

Images