Method for predicting drug-induced intrahepatic cholestasis based on multi-dimensional feature modeling

Abstract

The invention provides a method for predicting drug-induced intrahepatic cholestasis based on multi-dimensional feature modeling. The method comprises the following steps: collecting risk genes of intrahepatic cholestasis, and collecting intrahepatic cholestasis positive and negative drugs and targets thereof; establishing a disease cluster of intrahepatic cholestasis, and calculating a relative average shortest distance between the disease cluster and the drug target in the network to obtain network topology features; identifying a structural alarm of intrahepatic cholestasis and calculating Mordred descriptors of drug molecules to obtain molecular structure characteristics; carrying out modeling prediction by utilizing a machine learning algorithm; and finally, representing a domain of definition of the model by using a deep learning method. Compared with the traditional method, the method has the advantages that the big data artificial intelligence is integrated with the bioinformatics analysis technology, the drug characteristics are represented from the molecular structure and disease incidence of the drug, and a new method and a new thought are provided for predicting the complex toxicity of the drug.

Description

Technical Field

[0001] The present invention belongs to the field of pharmacy, and specifically relates to an artificial intelligence toxicity prediction method for drug-induced intrahepatic cholestasis based on network architecture analysis combined with chemical informatics. Background Art

[0002] Drug-induced liver injury (DILI) can be divided into cholestatic, hepatocellular and mixed types according to the type of tissue damage. Cholestatic DILI is the most common type of injury histologically, accounting for up to 50% of DILI cases, and can be divided into intrahepatic cholestasis and extrahepatic cholestasis according to the site of occurrence. Intrahepatic cholestasis is a liver disease caused by abnormalities in the structure and function of hepatocytes or capillary bile ducts due to various reasons, leading to abnormal bile formation, secretion and excretion. Cholestasis can lead to acute hepatotoxicity, increase the number of bile ducts and fibrosis, and then develop into cirrhosis, eventually leading to liver failure and the need for liver transplantation.

[0003] The pathogenesis of drug-induced intrahepatic cholestasis (DIIC) is complex, and there is currently a lack of effective therapeutic drugs. In addition, the systemic skin itching associated with intrahepatic cholestasis and long-term medication have caused great distress to patients. The prevalence of DIIC has restricted the launch of new drugs, and up to 70% of marketed drugs have been withdrawn by the FDA due to liver toxicity. Although drugs are increasingly becoming the cause of cholestasis, progress in accurately predicting and explaining mechanisms has not increased at the same rate. Traditional animal experiments to predict DIIC are not only time-consuming and costly, but also have species differences. Therefore, how to effectively and accurately predict DIIC remains an urgent problem to be solved in drug development and clinical drug evaluation. Summary of the invention

[0004] The purpose of the present invention is to provide a drug toxicity prediction method for quickly and accurately predicting DIIC. The present invention characterizes and models drug characteristics by considering the association between drug targets and structures and diseases through gene cluster proximity, structural alerts, and Mordred descriptors, and further uses deep learning to characterize the model definition domain, while analyzing the relationship between structural alerts and toxicity. The present invention can predict DIIC at low cost, high efficiency, and accurately, and provide a case reference for the prediction of other complex toxic endpoints.

[0005] The purpose of the present invention is achieved through the following technical solutions:

[0006] The method for predicting drug-induced intrahepatic cholestasis based on multidimensional feature modeling includes the following steps:

[0007] (1) Data collection and standardization

[0008] Risk genes for intrahepatic cholestasis were collected and standardized in the disease database using medical subject headings, and drugs with known DIIC labels and their targets were collected and standardized.

[0009] (2) Calculation and screening of drug characteristics

[0010] ① Clustering and functional annotation of DIIC risk genes to obtain disease clusters, and using the network proximity algorithm to calculate the relative average shortest path length between each disease cluster and the drug target in the PPI network, using Z dc Represents the topological characteristics of the biological network that characterizes the drug.

[0011] ② Use SARpy to identify the structural alerts of DIIC, calculate the IG value and PR value of each structural alert, and obtain the comprehensive value by weighting according to the IG value and PR value, which is expressed as SA. At the same time, use the Python package corresponding to Mordred to calculate the Mordred descriptor of each drug molecule and preprocess it. SA and Mordred descriptor jointly characterize the molecular structure characteristics of the drug.

[0012] ③ The biological network topological characteristics and molecular structure characteristics of the drugs were integrated, and recursive feature elimination cross validation (RFECV) was used to screen the feature set with the best accuracy under ten-fold cross validation.

[0013] (3) Model construction and evaluation

[0014] Multidimensional feature models were constructed using three machine learning algorithms, namely SVM, XGB and RF. For comparison, corresponding single-dimensional feature models were established using only Mordred descriptors. The accuracy (ACC), precision, sensitivity (SE), specificity (SP), F1 score (F1), Matthews correlation coefficient (MCC) and area under the receiver operating characteristic curve (AUC) of each model were calculated.

[0015] (4) Domain

[0016] KNN-GCN is used to predict the association between drugs and disease clusters. Drugs are assigned to different clusters based on the association to form a domain. Then, the drugs in each domain are modeled and the performance is evaluated separately. At the same time, random sampling with the same number of drugs in each cluster is used as a control group to model and evaluate the performance.

[0017] (5) Characteristics and mechanism of DIIC

[0018] The SHAP method was used to rank the importance of each feature and quantify the impact of the feature value on the positive and negative outputs of the model. Structural alerts that only existed in positive drugs in each DIIC cluster were screened, and the corresponding drugs were found. The targets of the drugs were used to overlap the network with the corresponding clusters, and finally the mechanism was explained through gene enrichment.

[0019] The medical subject terms referred to in step (1) of the present invention refer to the standardized terms for intrahepatic cholestasis in the MeSH database; disease databases include OMIM (https: / / www.omim.org / ), GeneCards (https: / / www.genecards.org / ), DisGeNET (https: / / disgenet.cn / ) and CTD (https: / / ctdbase.org / ) databases; drugs are collected in DILIrank and LiverTox and their labels are checked in metaADEDB (https: / / lmmd.ecust.edu.cn / online_services / metaadedb / ), and their targets are collected in DrugBank (https: / / go.drugbank.com / ) and DGIdb (https: / / www.dgidb.org / ); standardization refers to converting Uniprot IDs of all drug targets and disease risk genes into standard gene names in Uniprot (https: / / www.uniprot.org).

[0020] The functional annotation referred to in step (2) of the present invention refers to the gene ontology and pathway enrichment analysis of disease risk genes in the g:Profiler (https: / / biit.cs.ut.ee / gprofiler) database; the network proximity algorithm is used to calculate the relative average shortest distance between the drug target and the DIIC cluster, with Z dc It is calculated by formula 1-formula 2.

[0021]

[0022]

[0023] The IG of each structural alert represents information gain, which is calculated by Formula 3-Formula 5; PR represents the positive rate, which is calculated by Formula 6; aggregation is performed based on the IG value and PR value, and SA represents the aggregation value, which is calculated by Formula 7; the preprocessing of Mordred descriptors refers to deleting descriptors with most drugs as null values, variance less than 0.1, and Pearson coefficient greater than 0.9.

[0024] G(X)=-p1log2p1-p0log2p0 (Formula 3)

[0025]

[0026] IG(T)=H(X)-H(X|T) (Formula 5)

[0027]

[0028] SA=0.7PR+10IG (Formula 7)

[0029] The machine learning algorithm in step (3) of the present invention is called in Scikit-Learn, and the model performance is evaluated using the calculation results of Formula 8-Formula 15, where TP is true positive, TN is true negative, FP is false positive, FN is false negative, TPR is true positive rate, and TNR is true negative rate.

[0030] Accuracy = (TP + TN) / (TP + TN + FP + FN) (Formula 8)

[0031] P = Precision = TP / (TP + FP) (Formula 9)

[0032] Recall = SE = TP / (TP + FN) (Formula 10)

[0033] Specificity=SP=TN / (TN+FP) (Formula 11)

[0034] F1=(2×Precision×Recall) / (Precision+Recall) (Formula 12)

[0035]

[0036]

[0037] The KNN-GCN method described in step (4) of the present invention is to first generate edges between nodes by fitting the KNN model to construct a topological network structure, and then use the GCN algorithm to predict the association between disease clusters and drugs. GCN is used as a neural network layer, and the propagation definition from the 1st layer to the 1+1th layer is shown in Formula 15.

[0038]

[0039] The Shapley value in step (5) of the present invention is calculated by using the Shapley Additive Interpretation (SHAP) Python package to quantify the contribution of each input feature in the model; the attribution of each feature of the model is plotted into the SHAP bee colony, and the influence of the selected feature on the onset of DIIC hepatotoxicity is explored and explained based on the distribution of the high (or low) feature value of each feature in the coordinate quadrant; in addition, network overlap is performed by locating genes with ds≤2 in the PPI network, and gene enrichment to obtain key pathways is performed in the DAVID database (https: / / david.ncifcrf.gov / tools.jsp).

[0040] The present invention has the following beneficial effects and advantages:

[0041] 1. The present invention starts from disease genes and drug targets, taking into account the interconnections between disease genes and between disease genes and drug targets, rather than just considering these genes in isolation, and provides insights into the comprehensive nature of human pathogenesis.

[0042] 2. The present invention aims to improve the accuracy of hepatotoxicity prediction. Based on the latest molecular structure characteristics and the relationship characteristics of gene biological networks, it explores multi-dimensional feature modeling for predicting DIIC, a type of hepatotoxicity with a complex mechanism. Compared with traditional single-dimensional modeling, it has higher accuracy and more robust prediction performance.

[0043] 3. The present invention uses KNN to construct a network topology diagram of disease clusters, and GCN classifies drugs into different disease clusters related to DIIC for further exploration, extending from conventional toxicity modeling to mechanism modeling, which can improve the potential for prediction accuracy while helping to explore the intrinsic connection between drugs and disease onset.

[0044] 4. The method based on multidimensional feature modeling described in the present invention integrates multiple possible influencing factors in the drug pathogenesis, reveals the structure-drug-pathway-toxicity relationship, and provides case references for evaluating various complex toxicities of compounds. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 A schematic flow chart of the method of the present invention;

[0046] Figure 2 Clustering of disease genes for intrahepatic cholestasis;

[0047] Figure 3 Alerts on structures associated with intrahepatic cholestasis based on SARpy identification;

[0048] Figure 4 KNN results under five-fold cross validation;

[0049] Figure 5 Beeswarm plot of feature importance ranking based on SHAP. DETAILED DESCRIPTION

[0050] The following examples are only used to help understand the method of the present invention and its core concept. It should be pointed out that for ordinary technicians in this technical field, without departing from the principles of the present invention, several improvements and modifications can be made to the present invention, and these improvements and modifications also fall within the scope of protection of the claims of the present invention. The steps of the method for predicting drug-induced intrahepatic cholestasis based on multidimensional feature modeling in the embodiment of the present invention are as follows: Figure 1 shown.

[0051] Example 1

[0052] Prediction method of drug-induced intrahepatic cholestasis based on multidimensional feature modeling

[0053] 1. Methods and Steps

[0054] 1.1 Data Collection and Standardization

[0055] According to the MeSH medical subject headings "Intrahepatic cholestasis" and "Cholestasis, intrahepatic" for intrahepatic cholestasis as search keywords, disease genes were searched in OMIM (https: / / www.omim.org / ), GeneCards (https: / / www.genecards.org / ), DisGeNet (https: / / disgenet.cn / ), and CTD (https: / / ctdbase.org / ) databases, and Uniprot IDs were converted into standardized gene names in the Uniprot database (https: / / www.uniprot.org). Drugs in DILIrank and LiverTox were checked in the metaADEDB database (http: / / lmmd.ecust.edu.cn / metaadedb / ), and labels were defined as positive and negative for intrahepatic cholestasis, and drug targets were collected and standardized in DrugBank (https: / / go.drugbank.com / ) and DGIdb (https: / / www.dgidb.org / ).

[0056] 1.2 Calculation and screening of drug characteristics

[0057] ①Intrahepatic cholestasis risk genes were imported into the STRING database and the disease module was obtained using a confidence score > 0.9. The Louvain algorithm was then used to cluster the genes in the module. Functional enrichment of gene clusters was performed in the g:Profiler database (https: / / biit.cs.ut.ee / gprofiler). All genes in the complete DIIC module were used as the background gene set. Benjamin-Hochberg was used for multiple testing correction. After functional annotation of GO, KEGG, and REACTOME pathways, the DIIC disease clusters were obtained. The network proximity algorithm was used to calculate the relative average shortest distance (Z) between each disease cluster and the drug target. dc ), as the feature input of the model. dc Calculated by Formula 1-Formula 2.

[0058]

[0059]

[0060] ② Use SARpy software to extract drug structure alerts of DIIC, and use default parameters for extraction, that is, the number of atoms ranges from 2 to 18, and the minimum number of occurrences is 3. Then calculate the information gain value (IG) of each structure alert by formula 3-formula 5, calculate the positive rate value (PR) by formula 6, and finally calculate the SA aggregation value of each structure alert by formula 7. Then use the MordredPython package to calculate the 1D / 2D descriptors of all drug molecules, first delete the descriptors with null values ​​for most drugs, then delete the descriptors with a standard deviation less than 0.10, and finally delete the descriptors with a Pearson correlation coefficient greater than 0.90 to remove redundancy.

[0061] H(X)=-p1log2p1-p0log2p0 (Formula 3)

[0062]

[0063] IG(T)=H(X)-H(X|T) (Formula 5)

[0064]

[0065] SA=0.7PR+10IG (Formula 7)

[0066] ③ The multi-dimensional feature set is obtained by combining the above features and screening them using recursive feature elimination cross validation (RFECV). By obtaining the importance of each feature and recursively repeating this process on the pruned set, and then evaluating the accuracy of the selected feature set with 10-fold CV, the features with smaller contributions can be removed, and the feature set with the highest accuracy under 10-fold CV is selected as the DIIC multi-dimensional feature set. As a control, the single-dimensional feature set is the Mordred feature set with the highest accuracy under 10-fold CV selected in the same feature selection method.

[0067] 1.3 Model Construction and Evaluation

[0068] Three machine learning algorithms, extreme gradient boosting (XGB), random forest (RF) and support vector machine (SVM), were used in Scikit-Learn to establish a prediction model between DIIC hepatotoxicity and multidimensional features. The drugs collected by DILIrank were divided into training set and validation set at a ratio of 0.75:0.25, and the drugs collected by LiverTox were used as external test sets. The hyperparameters of the training set were optimized using a 10-fold CV grid search strategy to determine the best hyperparameter combination for each classifier. The multidimensional feature model and the single-dimensional feature model were modeled 100 times respectively. By performing 10-fold CV on the training set and applying the model to predict the validation set and external test set, the accuracy (ACC), precision (precision), recall (sensitivity, SE), specificity (SP), F1 score (F1), Matthews correlation coefficient (MCC) and area under the receiver operating characteristic curve (AUC) were calculated by formula 8-formula 14 to evaluate the model performance.

[0069] Accuracy = (TP + TN) / (TP + TN + FP + FN) (Formula 8)

[0070] Precision = TP / (Tp+FP) (Formula 9)

[0071] Recall = SE = TP / (TP + FN) (Formula 10)

[0072] Specificity=SP=TN / (TN+FP) (Formula 11)

[0073] F1=(2×Precision×Recall) / (Precision+Recall) (Formula 12)

[0074]

[0075]

[0076] 1.4 Domain

[0077] First, the structural similarity of drugs was calculated based on the Tanimoto similarity coefficient of ECFP. Second, the molecular function (MF) similarity of disease clusters was calculated using the R package GOSemSim. Drugs and diseases were represented as nodes, and the similarity matrix was represented as node features. Node labels were denoted by Z dc Value marker, negative Z dc When the value is 1 (association), positive Z dc The label is 0 (not associated) when the value is 0. Then, the edges between nodes are generated by fitting the KNN model to construct the topological network structure. Finally, the GCN algorithm is used to predict the association between disease onset clusters and drugs. As a neural network layer, the propagation definition of GCN from layer 1 to layer 1+1 is shown in Formula 15.

[0078]

[0079] Then, each disease cluster contains a different number of drugs, forming a set of related drugs within each cluster, representing the domain represented by KNN-GCN. The drugs belonging to DILIrank were randomly divided into training and validation sets at a ratio of 0.75:0.25, and the drugs belonging to LiverTox were used as external test sets. The DIIC cluster model was further constructed for the drugs included in each cluster. At the same time, the same number of positive drugs and negative drugs in each cluster were randomly selected from the original data set as the control group to establish the control model. The drugs in each disease cluster and each control group were randomly modeled 100 times to evaluate the KNN-GCN domain method.

[0080] 1.5 Characteristics and mechanism of DIIC

[0081] The Shapley value of each input feature was calculated by using the Shapley Additive Explanation (SHAP) Python package to quantify its contribution to the model output. The influence of the selected drug features on the incidence of DIIC hepatotoxicity was explored based on the distribution of high (or low) feature values ​​in the coordinate quadrants of the SHAP beeswarm plot. In addition, specific structural alerts that were only present in positive drugs were identified by analyzing the positive and negative occupancy of the selected SA among all drugs within each cluster. Subsequently, the network overlap of drug targets with the corresponding disease cluster was performed for the positive drugs containing this structural alert, and each drug target network was overlapped with the DIIC disease cluster module based on the same nodes in the PPI network. The ds between each specific risk gene of the DIIC disease cluster and its closest drug target in the biological network was calculated with a Bonferroni-corrected p-value threshold of 0.05, and then the genes with ds ≤ 2 were enriched to gain insight into the pathways through which structural alerts may affect DIIC disease.

[0082] 2. Results and Discussion

[0083] 2.1 Data Collection and Standardization

[0084] 456 drugs were collected in DILIrank and 295 drugs were obtained in LiverTox with known labels - DIIC hepatotoxicity (positive) or no DIIC hepatotoxicity (negative). Among them, 271 were DIIC positive and 480 were negative. In addition, 346 DIIC risk genes were obtained from OMIM, DisGeNET, GeneCards and CTD databases using the MeSH method.

[0085] 2.2 Calculation and screening of drug characteristics

[0086] ① First, the DIIC module was defined as 276 DIIC genes significantly located in the human interactome (Z = 12.6, p < 0.0001, STRING database), of which the largest connected component was composed of 250 risk genes (Z = 4.9, p = 1.4 × 10 -6 ). At the same time, when the modularity reached the maximum value of 0.63, the Louvain algorithm formed a total of 14 gene clusters. Only 8 statistically significant gene clusters with more than 10 genes were retained. Then the g:Profiler database was used to perform functional annotations on its biological pathways (p-value < 0.05), and a total of 8 DIIC clusters were obtained, such as Figure 2 The average shortest path length between each drug target and the DIIC cluster was calculated by the network proximity algorithm, and 8 Z dc , and use them to characterize the biological network topology characteristics between drugs and DIIC pathogenesis as the input of the model.

[0087] ② The 456 drugs selected from DILIrank were used as the training set, and the 19 structural alerts obtained using SARpy were as follows: Figure 3 As shown, the results show that the accuracy is 0.76, the sensitivity is 0.65, the specificity is 0.84, and the likelihood ratio (LR) of each structural alert is greater than 2, as shown in Table 1. The 295 drugs of LiverTox were used as the external test set to further detect these 19 structural alerts, and the prediction results were 0.69, 0.72, and 0.68, which verified that these 19 structural alerts were indeed closely related to DIIC. This means that for the collected drugs, their structural characteristics can be characterized by these 19 structural alerts that are highly correlated with DIIC. Then the IG and PR values ​​of each structural alert are calculated and weighted, and an aggregate value SA is assigned to each structural alert. The results are shown in Table 1. A matching search of the 19 structural alerts is performed on all drugs, and the SA value of the structural alert represents the SA characteristics of the drug containing the structural alert.

[0088] Table 1 19 structural alarms identified based on SARpy and their statistical results

[0089]

[0090] For Mordred molecular structure features, 1826 Mordred descriptors were calculated for 751 drugs. In order to remove redundant items, 1122 Mordred features remained after removing features with null values ​​for most drugs and standard deviations less than 0.10, and 422 Mordred features were retained after removing features with Pearson correlation coefficients greater than 0.90.

[0091] ③RFECV selected 8 multi-dimensional feature sets and 8 single-dimensional feature sets with the highest 10-fold CV accuracy. The multi-dimensional feature sets are SA, Zdc1, Zdc4, AATS4p, AATSC0i, GATS2d, GATS1v, SMR_VSA6, including 1 SA, 2 Z dc and 5 Mordred features. At the same time, for fair comparison, the selected single-dimensional features are 8 Mordred features: AATS1v, AATS4v, AATS4p, AATSC2s, GATS2d, GATS1v, GATS1i, and SLogP.

[0092] 2.3 Model Construction and Evaluation

[0093] Based on multi-dimensional and single-dimensional features, six models were established using three algorithms: RF, XGB and SVM. The parameter settings are detailed in Table 2.

[0094] Table 2 Hyperparameters of the models corresponding to single-dimensional features and multi-dimensional features

[0095]

[0096] Each model was randomly modeled 100 times, and the statistical results of its performance are shown in Table 3. In Table 3, the average ACC of the three models based on multidimensional features was always above 0.700, among which the RF model achieved an ACC of 0.801 in the validation set. However, the average ACC values ​​of the models based on single-dimensional features mostly remained below 0.700. Compared with the models based on single-dimensional features, the models based on multidimensional features showed superior predictive ability in all performance indicators in the three datasets. This obvious advantage emphasizes the enhanced ability of multidimensional features in capturing DIIC with complex disease bases.

[0097] Table 3 Performance results of multi-dimensional feature and single-dimensional feature models

[0098]

[0099] 2.4 Domain

[0100] The KNN-GCN method was used to classify 751 drugs into 8 DIIC clusters. The KNN results under 5-fold CV are shown in Figure 4 Although the performance is best when K = 1, this is usually meaningless. There is a rising peak when K = 9, which shows that when K = 9, KNN has better prediction ability than other values ​​of K. Therefore, K = 9 is selected to construct the topological network graph. The prediction results of GCN are shown in Table 4.

[0101] Table 4 5-fold CV results of GCN

[0102]

[0103] According to the KNN-GCN algorithm, 5502 drug cluster association pairs between 751 drugs and 8 DIIC incidence clusters were established. Since the prediction probability score exceeded 0.80 and the reliability was good, only 2671 pairs were used for the next model development analysis. At present, 751 drugs are divided into 8 cluster domains by the KNN-GCN classification method, and the specific number of drugs contained in each cluster is shown in Table 5. The RF model with the highest mean and the lowest standard deviation in 100 iterations with the best robustness and reliability is selected for cluster model construction. The results of 100 random modeling of 8 RF models between multidimensional features and DIIC hepatotoxicity based on 8 cluster definition domains are shown in Table 5. Among the 8 cluster models, the C3 model with ACC maintained above 0.800 has the best effect. Its ACC and AUC of the validation set are 0.903 and 0.837, respectively, and the ACC and AUC of the external test set are 0.890 and 0.810, respectively, indicating that the C3 model has strong robustness and predictive ability. In addition, according to the number of drugs in each domain, the same number of positive and negative drugs were selected as the corresponding control group in the domain. Then, the model was randomly built 100 times in the same way as the above cluster model. The results are listed in Table 5. It can be seen that the performance of each cluster model is significantly better than the control model.

[0104] Table 5 Performance results of cluster model and control model

[0105]

[0106]

[0107] *C1 represents lipid transport, C2 represents plasma remodeling, C3 represents metabolism of steroids, C4 represents drug metabolism-cytochrome P450, C5 represents inflammatory response, C6 represents apoptosis, C7 represents nucleosidephosphate biosynthetic process, C8 represents response to oxidative stress. n represents the number of the drugs in the cluster. a represents the control model, brepresents the cluster model.

[0108] 2.5 Characteristics and Mechanism of DIIC

[0109] The contribution of drug characteristics in the multidimensional feature RF model of 751 drugs is plotted into the SHAP bee colony separately, as shown in Figure 5 As shown in the SHAP value graph, SA features contribute the most, followed by Z dc 4 and Z dc 1, ranked second and fourth respectively, which shows that these two types of novel multi-dimensional features do have a significant effect on DIIC prediction. As can be seen from the figure, for SA, lower SHAP values ​​tend to be predicted as negative; on the contrary, higher SHAP values ​​are more likely to be predicted as positive, while Z dc Contrary to SA, it can be concluded that SA features are positively correlated with model output, Z dc 4 and Z dc 1 characteristic is negative correlation. It can be seen that the larger the SA value, the greater the DIIC toxicity risk. dc 4 and Z dc The lower the value of 1, the higher the risk of DIIC toxicity.

[0110] In addition, based on the fact that only structural alerts were contained in the positive drugs in each cluster, as shown in Table 1, drug-gene enrichment analysis was performed to explore the impact of the multidimensional characteristics of drugs on the onset of DIIC. SA3 (from Table 1, it can be seen that SA3 is furan), which may be closely related to DIIC. SA3 exists in the positive drug methoxsalen (CID: 4114). After its drug target overlapped with the C1 gene network, the enrichment results of ds≤2 genes showed that methoxsalen was related to transport regulation, among which ABCB11 (BSEP) was a key gene, which has been confirmed to be related to cholestasis. Methoxsalen may regulate the activity of ABCB11 through NR1H4 and NR0B2 in one step, affecting the normal transport of lipids, and NR1H4 and NR0B2 are neighbor genes of ABCB11. Previous studies have shown that methoxsalen is related to lipid transport by inhibiting ABCB11, causing bile acid accumulation and leading to hepatic cholestatic toxicity, which further confirmed our findings.

[0111] Similarly, SA7 was highlighted in C3 (which is a heteroatom chain containing nitrogen and sulfur as shown in Table 1), and SA7 may also be closely related to DIIC. Cimetidine (CID: 2756) is a positive drug containing SA7 in C3. The network overlap of cimetidine targets and genes in C3 showed that genes with ds≤2 were significantly enriched in the steroid biosynthesis pathway, which is a key step in steroid metabolism (C3). The enrichment results further showed that LSS, MSMO1 and SC5D regulate each other and participate in steroid biosynthesis together with hydroxymethylglutaryl coenzyme A (HMG-CoA). Previous studies have confirmed that cimetidine is indeed related to steroid biosynthesis, leading to DIIC.

[0112] In C4, SA3 and SA12 (shown as phenylthio in Table 1) that only appeared in positive drugs may have a strong correlation with DIIC. However, the drug-gene enrichment results of ash lactone (CID: 124039) containing SA3 and chlorpromazine (CID: 2726) containing SA12 did not show significant enrichment related to the C4 pathway. Therefore, their metabolites were further examined. Previous studies have shown that SA3 in ash lactone can be metabolized by CYP450 enzymes to cis-butene structure, leading to hepatotoxicity; SA12 in chlorpromazine can be metabolized by CYP450 enzymes to electrophilic intermediates, which then bind to glutathione, also leading to hepatotoxicity. Due to the structural differences between the parent drug and its metabolites when the drug is metabolized in the liver, it is often not the parent drug itself that produces hepatotoxicity, but its metabolites. This finding further emphasizes the importance of distinguishing the parent drug and its metabolites from hepatotoxicity for some drugs.

[0113] SA3 only appears in the positive drugs in C1 and C4, indicating that the SA3 structure may affect these two clusters and cause DIIC. dc 1 and Z dc 4 is negatively correlated with DIIC hepatotoxicity. By identifying drugs with significantly low values ​​of both Zdc1 and Zdc4 ​​and drugs with SA3 structures belonging to C1 or C4, 7 drugs: amiodarone (CID: 2157), methoxsalen (CID: 4114), dronedarone (CID: 208898), lapatinib (CID: 208908), nitrofurantoin (CID: 6604200), dantrolene (CID: 6914273), and cefuroxime (CID: 5479529) are indeed DIIC positive. Taking lapatinib (CID: 208908) as an example, the parent drug significantly inhibits ABCB1 and ABCG2 transporters, further affecting lipid excretion function, and has indeed been reported to show obvious DIIC. It is reported that its metabolite 8-hydroxylapatinib also interacts with active substances such as aldehydes or quinoneimines, and interacts with transporters, resulting in significant toxicity.

[0114] In summary, the use of multi-dimensional features to characterize drugs and the KNN-GCN method to classify the DIIC clusters to which each drug may belong can help improve the accuracy of DIIC hepatotoxicity prediction. This method can not only successfully predict the DIIC of compounds, but also logically infer the structure-pathway-toxicity relationship. These insights can serve as valuable methods for drug development and evaluation. The multi-dimensional feature model established in this patent is combined with the domain characterization to provide a new reference anchor pathway analysis model for exploring drug toxicity.

[0115] The above specific examples are only for understanding the technical solutions and core principles of the present invention, and are not intended to limit the scope of the present invention. For those skilled in the art, according to the core principles of the present invention, the conditions and parameters of the present examples that are changed and modified still belong to the protection scope of the present invention.

Claims

1. A method for predicting drug-induced intrahepatic cholestasis based on multidimensional feature modeling, comprising the following steps: (1) Data collection and standardization Risk genes for intrahepatic cholestasis were collected and standardized in the disease database using medical subject headings, and drugs with known DIIC labels and their targets were collected and standardized; (2) Calculation and screening of drug characteristics ① Clustering and functional annotation of DIIC risk genes to obtain disease clusters, and using the network proximity algorithm to calculate the relative average shortest path length between each disease cluster and the drug target in the PPI network, using Z dc It indicates the topological characteristics of biological networks that characterize drugs; ② Use SARpy to identify the structural alerts of DIIC, calculate the IG value and PR value of each structural alert, and obtain the comprehensive value by weighting the IG value and PR value, expressed as SA. At the same time, use the Python package corresponding to Mordred to calculate the Mordred descriptor of each drug molecule and preprocess it. SA and Mordred descriptors jointly characterize the molecular structure characteristics of the drug; ③ Comprehensively consider the biological network topological characteristics and molecular structure characteristics of the drugs, and use recursive feature elimination cross validation (RFECV) to screen the feature set with the best accuracy under ten-fold cross validation; (3) Model construction and evaluation Multidimensional feature models were constructed using three machine learning algorithms: SVM, XGB, and RF. For comparison, corresponding single-dimensional feature models were established using only the Mordred descriptor. The accuracy (ACC), precision (precision), recall (sensitivity, SE), specificity (SP), F1 score (F1), Matthews correlation coefficient (MCC), and area under the receiver operating characteristic curve (AUC) of each model were calculated. (4) Domain KNN-GCN is used to predict the association between drugs and disease clusters. Drugs are assigned to different clusters based on the association to form a domain. Then, the drugs in each domain are modeled and their performance is evaluated. At the same time, random sampling with the same number of drugs in each cluster is used as a control group to model and evaluate the performance. (5) Characteristics and mechanism of DIIC The SHAP method was used to rank the importance of each feature and quantify the impact of the feature value on the positive and negative outputs of the model. Structural alerts that only existed in positive drugs in each DIIC cluster were screened, and the corresponding drugs were found. The targets of the drugs were used to overlap the network with the corresponding clusters, and finally the mechanism was explained through gene enrichment.

2. The method according to claim 1, characterized in that The medical subject headings refer to the standardized terms for intrahepatic cholestasis in the MeSH database; disease databases include OMIM (https: / / www.omim.org / ), GeneCards (https: / / www.genecards.org / ), DisGeNET (https: / / disgenet.cn / ) and CTD (https: / / ctdbase.org / ) databases; drugs are collected in DILIrank and LiverTox and searched in metaADEDB Tags were checked in (https: / / lmmd.ecust.edu.cn / online_services / metaadedb / ), and their targets were collected in DrugBank (https: / / go.drugbank.com / ) and DGIdb (https: / / www.dgidb.org / ); Normalization refers to converting Uniprot IDs of all drug targets and disease risk genes into standardized gene names in Uniprot (https: / / www.uniprot.org).

3. The method according to claim 1, characterized in that Functional annotation refers to the gene ontology and pathway enrichment analysis of disease risk genes in the g:Profiler database (https: / / biit.cs.ut.ee / gprofiler); The network neighbor algorithm is used to calculate the relative average shortest distance between drug targets and DIIC clusters, with Z dc It is calculated by formula 1-formula 2.

4. The method according to claim 1, characterized in that The IG of the structural alert represents information gain, which is calculated by formula 3-formula 5; PR represents the positive rate, which is calculated by formula 6; aggregation is performed based on the IG value and the PR value, and SA represents the aggregation value, which is calculated by formula 7; the preprocessing of the Mordred descriptor refers to deleting most drugs as controls, descriptors with variances less than 0.1, and deleting descriptors with Pearson coefficients greater than 0.

9. H(X)=-P1log2p1-p0log2p0 (Formula 3) IG(T)=H(X)-H(X|T) (Formula 5) SA=0.7PR+10IG (Formula 7) 5. The method according to claim 1, characterized in that The machine learning algorithm is called in Scikit-Learn, and the model performance is evaluated using the calculation results of Formula 8-Formula 15, where TP is true positive, TN is true negative, FP is false positive, FN is false negative, TPR is true positive rate, and TNR is true negative rate. Accuracy = (TP + TN) / (TP + TN + FP + FN) (Formula 8) P = Precision = TP / (TP + FP) (Formula 9) Recall = SE = nP / (TP + FN) (Formula 10) Specificity=SP=TN / (TN+FP) (Formula 11) F1=(2×Precision×Recall) / (Precision+Recall) (Formula 12) 6. The method according to claim 1, characterized in that The KNN-GCN method first generates edges between nodes by fitting the KNN model to construct a topological network structure, and then uses the GCN algorithm to predict the association between disease clusters and drugs. GCN is used as a neural network layer, and the propagation definition from the 1st layer to the 1+1th layer is shown in Formula 15.

7. The method according to claim 1, characterized in that The Shapley values ​​were calculated using the Shapley Additive Explanations (SHAP) Python package to quantify the contribution of each input feature in the model; the attribution of each feature of the model was plotted into the SHAP bee colony, and the effect of the selected features on the onset of DIIC hepatotoxicity was explored and explained based on the distribution of high (or low) eigenvalues ​​of each feature in the coordinate quadrant; in addition, network overlap was performed by locating genes with ds≤2 in the PPI network, and gene enrichment to obtain key pathways was performed in the DAVID database (https: / / david.ncifcrf.gov / tools.jsp).

Images