A method of screening optimal descriptors and predicting classification of chemical activity

Abstract

The application discloses a method for screening optimal descriptors and predicting classification of chemical activity, comprising the following steps: (1) constructing a chemical data set; (2) calculating molecular descriptors; (3) preliminary screening of molecular descriptors; (4) dividing data samples; (5) screening molecular descriptors by using a combination of support vector machines and ant colony optimization algorithm; (6) counting the most frequent descriptors; (7) predicting the chemical activity by using support vector machines according to the most frequent descriptors; the application can make the classifier more robust by iteratively dividing the original training set and running the ant colony optimization algorithm; the initial setting of different path pheromones can make the ant colony optimization algorithm converge in a small number of iterations; the pheromone increment limit setting can avoid the algorithm falling into local optimum; and the statistical frequency can better quantify the importance of features and thus screen the optimal descriptor combination.

Description

Technical Field

[0001] This invention relates to cheminformatics, and more particularly to a method for screening optimal descriptors and predicting the activity classification of chemicals. Background Technology

[0002] Throughout their lives, humans are exposed to numerous different chemicals through various channels, including food, cleaning products, and pharmaceuticals. While some of these chemicals may be toxic, data on their potential harmful effects on human health remains insufficient. For thousands of pharmaceutical and environmental chemicals, there is a need to develop high-throughput methods to systematically assess their activity (toxicity). For example, the Tox21 project, a collaboration between the National Institute of Environmental Health Sciences (NIEHS), the National Toxicology Program (NTP), and the Environmental Protection Agency (EPA), aims to rapidly and effectively test whether certain chemicals may harm the human body and cause adverse health effects. It uses high-throughput screening (qHTS) to screen through a set of nuclear receptor (NR) and stress response (SR) pathways. The active molecules in Tox21 are those that can bind to one or more biochemical pathways and produce toxic effects in the human body. The toxic effects included in the Tox21 dataset refer to the effects of stress response (SR) and nuclear receptor (NR). Both of these effects are highly relevant to human health because activation of the SR pathway can lead to liver damage or cancer, while activation of the NR can disrupt the function of the endocrine system.

[0003] Quantitative structure-activity relationship (QSAR) models are an important tool for filling gaps in chemical environmental safety data. QSAR refers to the quantitative relationship between a molecule's microstructure and its physicochemical properties, environmental exposure behavior, in vivo toxicokinetics, and toxicological effect parameters (endpoints). In recent years, with the development of data analysis, machine learning, and artificial intelligence, techniques such as support vector machines and random forests have been used to build QSAR models. These models typically use molecular fingerprints or molecular descriptors to represent microstructure and physicochemical properties. Developing QSAR models for studying chemical toxicity has significantly reduced the cost of analyzing active (toxic) substances.

[0004] Currently, molecular descriptors are typically selected manually based on experimental properties. However, for high-dimensional molecules, this selection can easily lead to highly relevant or redundant descriptors, requiring significant expertise and causing inconvenience to the research. Therefore, manual selection has limitations and can result in unstable QSAR models. Furthermore, a large number of descriptors can reduce the general performance or predictive performance of the classifier. Therefore, selecting a small number of descriptors that truly affect the activity (toxicity) of chemicals is an important research direction in QSAR modeling. This can guide the structural optimization of existing active compounds and can also be used to determine the activity of new compound molecules. Summary of the Invention

[0005] Purpose of the invention: The purpose of this invention is to provide a method for screening the optimal descriptor and predicting the activity classification of chemicals to solve the problems existing in the above-mentioned technical background.

[0006] Technical solution: The method for screening optimal descriptors and predicting the classification of chemical activities according to the present invention includes the following steps:

[0007] (1) Constructing a chemical dataset: Constructing a chemical dataset using authoritative databases such as PubChem and ChEMBL, including: obtaining the SMILES code of the chemical, whether the chemical is active, and the molecular pathway endpoint.

[0008] (2) Calculate molecular descriptors: Based on the SMILES codes of the chemical activity dataset, the Mordred molecular descriptors of each chemical are calculated using the Mordred1.2 molecular descriptor calculator, including: composition descriptor, geometric descriptor, electrostatic topological descriptor, topological descriptor, and physicochemical descriptor.

[0009] (3) Preliminary screening of molecular descriptors: retain 2D descriptors and delete 3D descriptors, descriptors that cannot be calculated and descriptors with variance of 0.

[0010] (4) Data sample division: The chemicals in the dataset are randomly divided into 5 folds, with each 4 folds forming a training set and the other 1 fold forming a test set. This process is repeated 5 times.

[0011] (5) Use the ant colony optimization algorithm combined with support vector machine to screen molecular descriptors;

[0012] (6) The descriptor that appears most frequently in the statistical analysis;

[0013] (7) Predict chemical activity using support vector machine based on the most frequent descriptor.

[0014] Furthermore, step (5) involves dividing each 4-fold step (4) into a training set and a validation set within the ant colony optimization algorithm. Based on the set number of divisions, the training set and validation set are iteratively divided multiple times. Each division utilizes the ant colony optimization algorithm combined with a support vector machine to select the optimal descriptor combination for the current number of divisions. Specifically, this includes the following steps:

[0015] (51) Initialize parameters such as pheromone, number of ants, number of iterations, evaporation factor, and upper and lower bounds of pheromone increment;

[0016] (52) The first generation of ants selects descriptors based on probability;

[0017] (53) Evaluate using the fitness function based on the built-in support vectors;

[0018] (54) Retain the best combination of descriptors in the current algebra;

[0019] (55) Update pheromones;

[0020] (56) The next generation of ants repeats (53)-(55) based on the pheromone in step (55) and retains the optimal descriptor combination for the next generation;

[0021] (57) Compare the best descriptor combination of the current algebra with the best descriptor combination of the next generation until all algebras have run out, and select the best descriptor combination.

[0022] Furthermore, the probability calculation formula for step (52) is as follows:

[0023]

[0024] Where i represents the i-th station, i.e., the i-th molecular descriptor in the original training set, and j represents the path, which can be assigned a value of 0 or 1, indicating whether the corresponding descriptor is selected. τ ij This represents the pheromone concentration at station i on path j, while P ij Let and k represent the probability of being at the i-th station on the j-th path and the possible values ​​of path j, respectively.

[0025] Furthermore, the fitness function in step (53) is:

[0026] fitness=α*F-measure+β*G-mean+γ*MCCs.t.:α+β+γ=1 (2)

[0027] Furthermore, the formula for updating the pheromone in step (55) is:

[0028] τ ij (t+1)=ρ*τij (t)+Δτ ij (3)

[0029] Where p is the evaporation factor, referring to the proportion of pheromone reduction along each path, and Δτ ij This represents the increase in pheromone levels along some better paths; Δτ ij It can be obtained from the following formula:

[0030]

[0031] Where ant_n is the size of the ant population.

[0032] Furthermore, step (6) is as follows: the ant colony optimization algorithm in step (5) selects the best descriptor for each time and calculates the frequency, and selects the feature with the highest frequency, which is the final best descriptor combination.

[0033] Furthermore, step (7) involves: establishing a QSAR model based on the final optimal descriptor combination from step (6) using the support vector machine method, and predicting the chemical activity on the test set in step (4).

[0034] Beneficial Effects: Compared with existing technologies, this invention has the following significant advantages: The ant colony optimization algorithm adopts a direct-path structure, rather than an undirected graph model, which is more in line with the essence of feature selection; it makes the classifier more robust by iteratively dividing the original training set and running the ant colony optimization algorithm, and the setting of different path pheromone values ​​in the ant colony optimization algorithm can make the algorithm more inclined to select a smaller number of descriptors, thus achieving convergence quickly; it also sets upper and lower limits for pheromone increments to avoid the algorithm getting trapped in local optima; it divides the training set multiple times to count the frequency of the optimal descriptor, rather than the optimal descriptor combination at one time, which can quantify the importance of descriptors, better select the final optimal descriptor combination, and further help improve classification performance. Attached Figure Description

[0035] Figure 1 This is an overall flowchart of the screening descriptors and predicting chemical activity in this invention;

[0036] Figure 2 This is a framework diagram of the ant colony optimization algorithm for selecting descriptors in this invention;

[0037] Figure 3 The bar chart shows the frequency counts and rankings of the first 20 features in Embodiments 1 and 2 of this invention.

[0038] Figure 4 Line graphs comparing the results of selecting descriptors of different dimensions in Embodiments 1 and 2 of the present invention. Detailed Implementation

[0039] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0040] like Figures 1-2 As shown in the embodiment of the present invention, a method for screening optimal descriptors and predicting the classification of chemical activity includes the following steps:

[0041] (1) Construct a chemical dataset, specifically: construct a chemical dataset using authoritative databases such as PubChem and ChEMBL, including: obtaining the SMILES code of the chemical, whether the chemical is active, and the molecular pathway endpoint;

[0042] (2) Calculate molecular descriptors, specifically: based on the chemical activity dataset and SMILES code, calculate the Mordred molecular descriptors for each chemical using Mordred1.2 molecular descriptors, including: composition descriptor, geometric descriptor, electrostatic topological descriptor, topological descriptor, and physicochemical descriptor.

[0043] (3) Preliminary screening of molecular descriptors: retain 2D descriptors and delete 3D descriptors, as well as descriptors that cannot be calculated and have a variance of 0.

[0044] (4) The data sample is divided as follows: the chemicals in the dataset are randomly divided into 5 folds, each 4 folds is a training set, and the other 1 fold is a test set. This process is repeated 5 times.

[0045] (5) Using the ant colony optimization algorithm to screen molecular descriptors, specifically: dividing each 4-fold step in step (4) into a training set and a validation set within the ant colony optimization algorithm, and iterating the training set and validation set multiple times according to the set number of divisions. Each division can use the ant colony optimization algorithm to select the optimal combination of descriptors for the current number of divisions, specifically including the following steps:

[0046] (51) Parameters such as path initialization pheromone, number of ants, number of iterations, evaporation factor, and upper and lower bounds of pheromone increment are as follows: During the descriptor selection process, there are paths 0 and 1. We set the pheromone for path 0 to be greater than that for path 1. This makes the selection of descriptors more inclined to select a smaller number of descriptors. The ant colony optimization algorithm can achieve convergence without setting too many iterations. At the same time, by setting the upper and lower bounds of pheromone concentration, we can avoid the algorithm from getting trapped in local optima.

[0047] (52) The first generation of ants selects descriptors based on probability, and the probability calculation formula is:

[0048]

[0049] Where i represents the i-th station, i.e., the i-th molecular descriptor in the original training set, and j represents the path, which can be assigned a value of 0 or 1, indicating whether the corresponding descriptor is selected. τ ij This represents the pheromone concentration at station i on path j, while P ij Let and k represent the probability of the i-th station being on the j-th path and the possible values ​​of path j, respectively;

[0050] (53) Evaluate using the built-in support vectors and the fitness function, which is:

[0051] fitness=α*F-measure+β*G-mean+γ*MCC st:α+β+γ=1 (2)

[0052] (54) Retain the best combination of descriptors in the current algebra;

[0053] (55) Update the pheromone, wherein the formula for updating the pheromone is:

[0054] τ ij (t+1)=ρ*τ ij (t)+Δτ ij (3)

[0055] Where p is the evaporation factor, referring to the proportion of pheromone reduction along each path, and Δτ ij This represents the increase in pheromone levels along some better paths; Δτ ij It can be obtained from the following formula:

[0056]

[0057] Where ant_n is the size of the ant population;

[0058] (56) The next generation of ants repeats (53)-(55) based on the pheromone in step (55) and retains the optimal descriptor combination for the next generation;

[0059] (57) Compare the best descriptor combination of the current algebra with the best descriptor combination of the next generation until all algebras have run out, and select the best descriptor combination.

[0060] (6) Statistically determine the descriptor that appears most frequently. Specifically, in step (5), the ant colony optimization algorithm selects the optimal descriptor for each step and calculates the frequency. The descriptor with the highest frequency is selected, which is the final optimal descriptor combination.

[0061] (7) Based on the descriptor with the highest frequency, the chemical activity is predicted using a support vector machine. Specifically, based on the final optimal combination of descriptors from step (6) and using the support vector machine method, a QSAR model is established, and the chemical activity is predicted on the test set in step (4).

[0062] Example 1

[0063] Given a Tox21 dataset, the molecular pathway endpoint is the androgen receptor MDA-kb2 cell line (NR-AR), containing 7044 molecules (represented by SMILES codes), of which 6743 are inactive molecules and 301 are active molecules, with a class imbalance rate of 4.3%, as shown in the table below:

[0064]

[0065] The Mordred molecular descriptors for each molecule are calculated using the SMILES codes of the chemicals and the Mordred 1.2 molecular descriptor calculator, totaling 1613 descriptors.

[0066] Descriptors that could not be computed and had a variance of 0 were removed, resulting in 672 descriptors.

[0067] The 7044 molecules are divided into training and test sets. This involves randomly partitioning the dataset into five folds, with each fold comprising a training set and the remaining fold a test set, repeating this process five times. In the ant colony optimization algorithm, each fold of the training set is further divided into training and validation sets. The best feature subset is selected from the training set, and the best feature subset is evaluated using a fitness function based on the support vector machine (SVM) in the validation set. The initialization parameters for the ant colony optimization algorithm with SVM are set as shown in the table below.

[0068]

[0069]

[0070] This study uses an ant colony optimization algorithm to screen the 672 initially selected descriptors. The algorithm uses 30 ants, 10 generations, and divides the training set 40 times. During descriptor selection, there are paths 0 and 1. We set the pheromone concentration for path 0 to be higher than that for path 1. This biases the selection towards fewer descriptors, allowing the ant colony optimization algorithm to converge without requiring too many iterations. Furthermore, we set upper and lower bounds for the pheromone concentration to prevent the algorithm from getting trapped in local optima. Finally, we statistically analyze the frequency of descriptor occurrences. Figure 3 As shown in a, the 20 most frequently occurring descriptors are filtered out as shown in the table below:

[0071]

[0072] To validate the selected descriptors for the ant colony algorithm, a support vector machine was used on the test set for verification. The results were compared only with the 672 descriptors selected in the initial screening. The results are as follows:

[0073] Number of descriptors F-measure G-mean MCC 20 0.603 0.687 0.617 672 0.552 0.647 0.573

[0074] As can be seen from the table above, the features selected by the ant colony algorithm, when evaluated by our evaluation function, have significantly higher metrics than the support vector machine QSAR model after the initial screening.

[0075] To make the features selected by the ant colony optimization algorithm more meaningful, we controlled for variables and compared the top 5, 10, 20, 30, 50, and 100 molecule descriptors in the same run. The results are as follows:

[0076] Number of descriptors F-measure G-mean MCC 5 0.574 0.670 0.587 10 0.597 0.686 0.608 20 0.603 0.687 0.617 30 0.604 0.688 0.618 50 0.609 0.694 0.621 100 0.608 0.694 0.621

[0077] Based on the results, it was found that for cell lines whose molecular pathway endpoint is the androgen receptor MDA-kb2, such as Figure 4 As shown in a, b, and c, the F-means, G-means, and MCC values ​​are highest at 50 features. However, the metrics at 10, 20, and 30 features are very close to the evaluation metrics at 50 features. Considering that the first 50 features include the calculated features at 10, 20, and 30, it can be concluded that the optimal descriptors in the 10-50 range selected by the ant colony optimization algorithm can be used to build the QSAR model.

[0078] Example 2

[0079] Given a Tox21 dataset, the molecular pathway endpoint is the androgen receptor GeneBLAzer AR-UAS-bla-GripTite cell line (NR-AR-LBD), containing 6572 molecules (represented by SMILES codes), of which 6349 are inactive molecules and 223 are active molecules, with a class imbalance rate of 3.4%, as shown in the table below:

[0080]

[0081] The Mordred molecular descriptors for each molecule are calculated using the SMILES codes of the chemicals and the Mordred 1.2 molecular descriptor calculator, totaling 1613 descriptors.

[0082] Descriptors that could not be computed and had a variance of 0 were removed, resulting in 669 descriptors.

[0083] The 6572 molecules were divided into training and test sets. This involved randomly partitioning the dataset into five folds, with each fold comprising a training set and the remaining fold a test set, repeating this process five times. In the ant colony optimization algorithm, each fold of the training set was further divided into training and validation sets. The best feature subset was selected from the training set, and the best feature subset was evaluated using a fitness function based on the support vector machine (SVM) in the validation set. The initialization parameters for the ant colony optimization algorithm with SVM are shown in the table below.

[0084]

[0085] This study uses an ant colony optimization algorithm to screen the 669 descriptors initially selected. The algorithm uses 30 ants, 10 generations of iterations, and divides the training set 40 times. During descriptor selection, there are paths 0 and 1. We set the pheromone concentration for path 0 to be higher than that for path 1. This biases the selection of fewer descriptors, allowing the ant colony optimization algorithm to converge without requiring too many iterations. Furthermore, we set upper and lower bounds for the pheromone concentration to prevent the algorithm from getting trapped in local optima. Finally, we statistically analyze the frequency of descriptor occurrences. Figure 3 As shown in b, the 20 most frequently occurring descriptors are filtered out as shown in the table below:

[0086]

[0087] To validate the selected descriptors for the ant colony algorithm, a support vector machine was used for validation on the test set. The results were compared only with the 669 descriptors selected in the initial screening. The results are as follows:

[0088] Number of descriptors F-measure G-mean MCC 20 0.617 0.729 0.614 669 0.573 0.679 0.582

[0089] As can be seen from the table above, the features selected by the ant colony algorithm, when evaluated by our evaluation function, have significantly higher metrics than the support vector machine QSAR model after the initial screening.

[0090] To make the features selected by the ant colony optimization algorithm more meaningful, we controlled for variables and compared the top 5, 10, 20, 30, 50, and 100 molecule descriptors in the same run. The results are as follows:

[0091] Number of descriptors F-measure G-mean MCC 5 0.572 0.718 0.564 10 0.616 0.729 0.613 20 0.617 0.729 0.614 30 0.608 0.720 0.608 50 0.607 0.720 0.606 100 0.613 0.721 0.613

[0092] Based on the results, it was found that for the molecular pathway endpoint androgen receptor GeneBLAzer AR-UAS-bla-GripTite cell line dataset, such as Figure 4As shown in d, e, and f, the F-measure, G-mean, and MCC are at the 20-feature range, which means that the ant colony optimization algorithm selects 20 features, producing the best results. However, the values ​​at the 10, 30, and 50-feature ranges are also very close to the optimal solution. Therefore, it can be concluded that the optimal descriptors in the 10-50 range selected by the ant colony optimization algorithm can be used to build the QSAR model.

Claims

1. A method for screening optimal descriptors and predicting the classification of chemical activity, characterized in that, Includes the following steps: (1) Constructing a chemical dataset: Using authoritative databases such as PubChem and ChEMBL, a chemical dataset is constructed, including: obtaining the SMILES codes of chemicals, whether chemicals are active, and molecular pathway endpoints; (2) Calculate molecular descriptors: Based on the SMILES codes of the chemical dataset, use the Mordred1.2 molecular descriptor calculator to calculate the Mordred molecular descriptors of each chemical, including: composition descriptor, geometric descriptor, electrostatic topological descriptor, topological descriptor, and physicochemical descriptor; (3) Preliminary screening of molecular descriptors: retain 2D descriptors and delete 3D descriptors, descriptors that cannot be calculated and descriptors with variance of 0; (4) Data sample division: The chemicals in the dataset are randomly divided into 5 folds, with each 4 folds forming a training set and the other 1 fold forming a test set. This process is repeated 5 times. (5) Use the ant colony optimization algorithm with support vector machine to screen molecular descriptors; divide each 4-fold step in step (4) into a training set and a validation set within the ant colony optimization algorithm. According to the set number of divisions, iterate the division of the training set and the validation set multiple times. Each division can use the ant colony optimization algorithm with support vector machine to select the optimal combination of descriptors for the current number of divisions. Specifically, the following steps are included: (51) Initialize path pheromone, number of ants, number of iterations, evaporation factor, and upper and lower bounds of pheromone increment parameters; (52) The first generation of ants selects descriptors based on probability; (53) Evaluate using the fitness function based on the built-in support vectors; the fitness function is: (2)； (54) Retain the best combination of descriptors in the current algebra; (55) Update pheromones; the formula for updating pheromones is: (3)； in, The evaporation factor refers to the percentage of pheromone reduction along each pathway. This represents the increase in pheromones along some better paths; Obtained from the following formula: (4)； in, The size of the ant population; (56) The next generation of ants repeats (53)-(55) based on the pheromone in step (55) and retains the optimal descriptor combination for the next generation; (57) Compare the best descriptor combination of the current algebra with the best descriptor combination of the next generation until all algebras have run out, and select the best descriptor combination. (6) In the statistical step (5), the ant colony optimization algorithm selects the best descriptor for each time and calculates the frequency. The feature with the highest frequency is selected, which is the final best descriptor combination. (7) Based on the final optimal descriptor combination in step (6), a QSAR model is established using the support vector machine method, and the chemical activity is predicted on the test set in step (4).

2. The method for screening optimal descriptors and predicting the activity classification of chemicals according to claim 1, characterized in that, In step (52), the probability calculation formula is: (1)； in Representing the The site, i.e., the site in the original training set. A molecular descriptor, The path is represented by a value of 0 or 1, indicating whether the corresponding descriptor is selected. Then it represents the first The site at the The pheromone concentration along the path, and and They represent the first time in the second month. The site at the Probability and path on a path The value of .

Images