Method for predicting intrahepatic cholestasis of drugs based on multi-dimensional feature modeling

Abstract

The application provides a method for predicting intrahepatic cholestasis of drugs based on multi-dimensional feature modeling. The application collects risk genes of intrahepatic cholestasis, collects positive and negative drugs of intrahepatic cholestasis and their targets, establishes a disease onset cluster of intrahepatic cholestasis, calculates the relative average shortest distance of the disease onset cluster and the drug target in the network to obtain network topology features, identifies the structural alarm of intrahepatic cholestasis and calculates the Mordred descriptor of the drug molecule to obtain molecular structure features, uses a machine learning algorithm for modeling and prediction, and finally uses a deep learning method to represent the definition domain of the model. Compared with the traditional method, the application integrates bioinformatics analysis technology by using big data artificial intelligence, characterizes the drug from the molecular structure and disease onset, provides a new method and new idea for the complex toxicity prediction of drugs.

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 based on the type of tissue damage. Cholestatic DILI is the most common histological type of injury, accounting for up to 50% of DILI cases, and can be divided into intrahepatic and extrahepatic cholestasis based on the site of occurrence. Intrahepatic cholestasis is a liver disease caused by abnormalities in the structure and function of hepatocellular or bile capillary structures 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, ultimately causing 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. Furthermore, the systemic itching associated with DIIC and the long-term use of medications cause significant distress to patients. The prevalence of DIIC has restricted the availability of new drugs, and up to 70% of marketed drugs have been withdrawn by the FDA due to hepatotoxicity. Although drugs are increasingly becoming a cause of cholestasis, progress in accurately predicting and explaining the mechanisms has not increased at the same rate. Traditional animal experiments to predict DIIC are not only time-consuming and costly, but also subject to species differences. Therefore, how to effectively and accurately predict DIIC remains an urgent issue in drug development and clinical drug evaluation. Summary of the Invention

[0004] The present invention aims to provide a drug toxicity prediction method that rapidly and accurately predicts DIIC. This method characterizes and models drug features using gene cluster proximity, structural alerts, and Mordred descriptors, comprehensively considering drug targets and the relationship between structure and disease. Deep learning is further used to characterize the model domain and analyze the relationship between structural alerts and toxicity. This method enables cost-effective, efficient, and accurate prediction of DIIC, providing a reference for the prediction of other complex toxicity 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 using medical subject headings in the disease database. Drugs with known DIIC signatures and their targets were also 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] ② SARpy was used to identify the structural alerts of DIIC, and the IG value and PR value of each structural alert were calculated. The comprehensive value was obtained by weighting the IG value and PR value, which was expressed as SA. At the same time, the Mordred descriptor of each drug molecule was calculated and preprocessed using the Python package corresponding to Mordred. SA and Mordred descriptor jointly characterized the molecular structural characteristics of the drug.

[0012] ③ 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.

[0013] (3) Model construction and evaluation

[0014] 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.

[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, thus forming a domain. The drugs within each domain are then modeled and evaluated for performance. A random sample of drugs with the same number of drugs within each cluster is also used as a control group for modeling and performance evaluation.

[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 headings referred to in step (1) of the present invention refer to standardized terms for intrahepatic cholestasis in the MeSH database; disease databases include OMIM, GeneCards, DisGeNET and CTD databases for retrieval; drugs are collected in DILIrank and LiverTox and their tags are checked in metaADEDB, and their targets are collected in DrugBank and DGIdb; standardization refers to converting the Uniprot IDs of all drug targets and disease risk genes into standard gene names in Uniprot.

[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 database; the network proximity algorithm is used to calculate the relative average shortest distance between drug targets and DIIC disease clusters, with Z dc Indicates that 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 the aggregation value is represented by SA, 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] H(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 to 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] 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 layer 1 to layer 1+1 is shown in Formula 15.

[0037]

[0038] 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 features on the pathogenesis of DIIC hepatotoxicity is explored and explained based on the distribution of high (or low) eigenvalues ​​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.

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

[0040] 1. The present invention starts from disease genes and drug targets, considers the connections 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.

[0041] 2. This 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.

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

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

[0044] Figure 1 Schematic diagram of the method flow of the present invention;

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

[0046] Figure 3 Alerts based on structures associated with intrahepatic cholestasis identified by SARpy;

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

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

[0049] 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.

[0050] Example 1

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

[0052] 1. Methods and Steps

[0053] 1.1 Data Collection and Standardization

[0054] We searched for disease genes in the OMIM, GeneCards, DisGeNet, and CTD databases using the MeSH terms "Intrahepatic cholestasis" and "Cholestasis, intrahepatic" for intrahepatic cholestasis. Uniprot IDs were converted to standardized gene names in the Uniprot database. We also searched for drugs in the metaADEDB databases (DILIrank and LiverTox), defined positive and negative labels for intrahepatic cholestasis, and collected and standardized drug targets in DrugBank and DGIdb.

[0055] 1.2 Calculation and screening of drug characteristics

[0056] ① Intrahepatic cholestasis risk genes were imported into the STRING database and a confidence score > 0.9 was used to obtain disease modules. The Louvain algorithm was then used to cluster the genes in the modules. Functional enrichment of gene clusters was performed in the g:Profiler database. All genes in the complete DIIC module were used as the background gene set. Benjamin-Hochberg correction for multiple testing was performed. After functional annotation of GO, KEGG, and REACTOME pathways, 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.

[0057]

[0058] ② SARpy software was used to extract drug structure alerts for DIIC using default parameters, i.e., an atom number range of 2 to 18 and a minimum number of occurrences of 3. The information gain (IG) value for each structure alert was then calculated using Equations 3-5, the positive rate (PR) value was calculated using Equation 6, and finally, the SA aggregate value for each structure alert was calculated using Equation 7. The Mordred Python package was then used to calculate 1D / 2D descriptors for all drug molecules. Descriptors with a majority of null values ​​were first removed, followed by descriptors with a standard deviation less than 0.10, and finally by descriptors with a Pearson correlation coefficient greater than 0.90 to eliminate redundancy.

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

[0060]

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

[0062]

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

[0064] ③ The multidimensional 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, the accuracy of the selected feature set is evaluated using 10-fold CV. Features with less contribution are removed, and the feature set with the highest accuracy under 10-fold CV is selected as the DIIC multidimensional feature set. As a control, the single-dimensional feature set is the Mordred feature set with the highest accuracy under 10-fold CV, selected using the same feature selection method.

[0065] 1.3 Model Construction and Evaluation

[0066] Three machine learning algorithms, extreme gradient boosting (XGB), random forest (RF), and support vector machine (SVM), were used in Scikit-Learn to establish a predictive model linking DIIC hepatotoxicity with multidimensional features. The drugs collected by DILIrank were divided into training and validation sets at a ratio of 0.75:0.25, and the drugs collected by LiverTox served as an external test set. Hyperparameter optimization was performed on the training set using a 10-fold CV grid search strategy to determine the optimal hyperparameter combination for each classifier. Multidimensional and unidimensional feature models were each built 100 times. Ten-fold CV was performed on the training set and the resulting models were used to predict the validation and external test sets. Model performance was evaluated by calculating accuracy (ACC), precision (SE), recall (SE), specificity (SP), F1 score (F1), Matthews correlation coefficient (MCC), and area under the receiver operating characteristic curve (AUC) using Equations 8-14.

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

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

[0069] Recall = SE = TP / (TP + FN) (Formula 10) Specificity = SP = TN / (TN + FP) (Formula 11)

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

[0071]

[0072] 1.4 Domain

[0073] 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 represented by Z dc Value marker, negative Z dc When the value is 1 (association), Z dc When the value is 0, the label is 0 (no association). Then, by fitting the KNN model, edges between nodes are generated 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 of the GCN from layer 1 to layer 1+1 is defined as shown in Equation 15.

[0074]

[0075] Each disease cluster then contained a varying number of drugs, forming a set of related drugs within each cluster, representing the domain represented by the KNN-GCN. Drugs belonging to DILIrank were randomly divided into training and validation sets at a ratio of 0.75:0.25. Drugs belonging to LiverTox were used as an external test set. A DIIC cluster model was further constructed for the drugs included in each cluster. Simultaneously, the same number of positive and negative drugs from each cluster were randomly sampled from the original dataset as control groups to establish control models. Drugs in each disease cluster and each control group were randomly modeled 100 times to evaluate the KNN-GCN domain approach.

[0076] 1.5 Characteristics and Mechanism of DIIC

[0077] The Shapley Additive Explanation (SHAP) Python package was used to calculate the Shapley value of each input feature to quantify its contribution to the model output. The impact of 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. Furthermore, specific structural alerts present only in positive drugs were identified by analyzing the positive and negative occupancy of selected SAs across all drugs within each cluster. Subsequently, the drug-target network of positive drugs containing these structural alerts was overlapped with the corresponding disease cluster. Each drug-target network was then overlapped with the DIIC disease cluster module based on the same nodes in the PPI network. The ds between each specific risk gene in 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. Genes with ds ≤ 2 were then enriched to gain insight into pathways where structural alerts may influence DIIC disease.

[0078] 2. Results and Discussion

[0079] 2.1 Data Collection and Standardization

[0080] 456 drugs were collected in DILIrank, and 295 drugs were identified in LiverTox with known labels—either DIIC-positive (positive) or non-DIIC-negative (negative). Of these, 271 were DIIC-positive and 480 were DIIC-negative. Furthermore, 346 DIIC-risk genes were identified using the MeSH approach from the OMIM, DisGeNET, GeneCards, and CTD databases.

[0081] 2.2 Calculation and screening of drug characteristics

[0082] ① First, the DIIC module was defined as 276 DIIC genes significantly mapped to 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 ). Meanwhile, when the modularity reached its maximum value of 0.63, the Louvain algorithm generated a total of 14 gene clusters. Only eight statistically significant gene clusters with more than 10 genes were retained. The g:Profiler database was then used to perform functional annotation of their biological pathways (p-value < 0.05), resulting in a total of eight DIIC clusters, 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.

[0083] ② 456 drugs selected from DILIrank were used as training sets, and 19 structural alerts obtained using SARpy were used. Figure 3 The results showed an accuracy of 0.76, a sensitivity of 0.65, and a specificity of 0.84, with the likelihood ratio (LR) for each structure alert greater than 2, as shown in Table 1. These 19 structure alerts were further tested using the 295 drugs in the LiverTox external test set, resulting in predictions with an accuracy of 0.69, a sensitivity of 0.72, and a specificity of 0.68, confirming that these 19 structure alerts are indeed closely associated with DIIC. This means that the structural characteristics of the collected drugs can be characterized by these 19 structure alerts that are highly correlated with DIIC. The IG and PR values ​​for each structure alert were then calculated and weighted, and an aggregate value SA was assigned to each structure alert. The results are shown in Table 1. A match search for the 19 structure alerts was performed on all drugs, and the SA value of each structure alert represented the SA characteristics of the drug containing the structure alert.

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

[0085]

[0086] For Mordred molecular structure features, 1826 Mordred descriptors were calculated for 751 drugs. To remove redundant terms, features with null values ​​for most drugs and standard deviations less than 0.10 were removed, resulting in 1122 remaining Mordred features. After removing features with Pearson correlation coefficients greater than 0.90, 422 Mordred features were retained.

[0087] ③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 single-dimensional features selected are 8 Mordred features: AATS1v, AATS4v, AATS4p, AATSC2s, GATS2d, GATS1v, GATS1i, and SLogP.

[0088] 2.3 Model Construction and Evaluation

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

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

[0091]

[0092] Each model was randomly tested 100 times, and the statistical results of their performance are shown in Table 3. In Table 3, the average ACC of the three models based on multidimensional features was consistently above 0.700, with the RF model achieving 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 demonstrated superior predictive ability across all performance metrics in the three datasets. This clear advantage emphasizes the enhanced ability of multidimensional features in capturing DIIC with complex disease bases.

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

[0094]

[0095] 2.4 Domain

[0096] 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 K = 1 shows the best performance, this is generally meaningless. A rising peak appears when K = 9, indicating that KNN has better prediction ability when K = 9 than other values ​​of K. Therefore, K = 9 was selected to construct the topological network graph. The prediction results of GCN are shown in Table 4.

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

[0098]

[0099] Using the KNN-GCN algorithm, 5502 drug cluster association pairs were established between 751 drugs and eight DIIC incidence clusters. Because the predicted probability scores exceeded 0.80, indicating good reliability, only 2671 pairs were used for the next round of model development and analysis. The KNN-GCN classification method now categorizes the 751 drugs into eight cluster domains. The specific number of drugs in each cluster is shown in Table 5. The RF model with the highest mean and lowest standard deviation across 100 iterations, demonstrating the best robustness and reliability, was selected for cluster model construction. The results of 100 random iterations of eight RF models linking multidimensional features to DIIC hepatotoxicity based on the eight cluster domains are shown in Table 5. Among the eight cluster models, the C3 model, which maintained an ACC above 0.800, performed best, with an ACC and AUC of 0.903 and 0.837 on the validation set and 0.890 and 0.810 on the external test set, respectively, demonstrating the C3 model's robustness and predictive power. In addition, based on the number of drugs in each domain, we selected the same number of positive and negative drugs as the corresponding control group in that domain. We then randomly modeled 100 times in the same manner as the cluster model above. The results are listed in Table 5, showing that the performance of each cluster model is significantly better than the control model.

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

[0101]

[0102]

[0103] *C1 represents regulation oflipid metabolism by PPARalpha,C2represents plasma lipoprotein remodeling,C3represents metabolism of steroids,C4 represents drugmetabolism-cytochrome P450,C5 represents inflammatory response,C6 represents apoptosis,C7 represents nucleosidephosphatebiosynthetic process,C8represents response to oxidativestress.n represents the number of the drugs in the cluster. a represents the control model,b represents the cluster model.

[0104] 2.5 Characteristics and Mechanism of DIIC

[0105] The contributions of drug characteristics in the multidimensional feature RF model of 751 drugs are plotted into the SHAP bee colony, as shown in Figure 5 As shown in the SHAP value graph, SA features have the largest contribution, 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, and Z dc In contrast 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.

[0106] Furthermore, as shown in Table 1, drug-gene enrichment analysis was performed to explore the impact of multidimensional drug features on the pathogenesis of DIIC, given that positive drugs within each cluster contained only structural alerts. SA3 (as shown in Table 1, SA3 is a furan) is associated with C1 and may be closely related to DIIC. SA3 is present in the positive drug methoxsalen (CID: 4114). After overlapping its drug targets with the C1 gene network, enrichment results for genes with a ds ≤ 2 revealed that methoxsalen is involved in transport regulation, with PPARα being a key gene, which has been linked to cholestasis. Methoxsalen may regulate lipid transport through PPARα, thereby affecting normal lipid metabolism. Previous studies have shown that methoxsalen, by inhibiting PPARα, causes bile acid accumulation, leading to hepatic cholestatic toxicity, further corroborating our findings.

[0107] Similarly, SA7 (a nitrogen- and sulfur-containing heteroatom chain, as shown in Table 1) is prominently identified in C3 and may also be closely associated with DIIC. Cimetidine (CID: 2756) is a positive drug in C3 for SA7. Network overlap between cimetidine targets and genes in C3 revealed that genes with a ds ≤ 2 were significantly enriched in the steroid biosynthesis pathway, a key step in steroid metabolism (C3). The enrichment results further suggest 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 involved in steroid biosynthesis, leading to DIIC.

[0108] In C4, SA3 and SA12 (shown as phenylthio groups in Table 1), which appear only in positive drugs, may have a strong association with DIIC. However, drug-gene enrichment results for salvinolactone (CID: 124039), which contains SA3, and chlorpromazine (CID: 2726), which contains SA12, did not show significant enrichment related to the C4 pathway. Therefore, their metabolites were further examined. Previous studies have shown that SA3 in salvinolactone can be metabolized by CYP450 enzymes to a cis-butene structure, leading to hepatotoxicity; SA12 in chlorpromazine can be metabolized by CYP450 enzymes to an electrophilic intermediate, which then conjugates to glutathione, also leading to hepatotoxicity. Due to the structural differences between the parent drug and its metabolites during hepatic metabolism, hepatotoxicity is often not caused by the parent drug itself, but by its metabolites. This finding further emphasizes the importance of distinguishing between the parent drug and its metabolites in terms of hepatotoxicity for certain drugs.

[0109] SA3 only appeared 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 Zdc1 and Zdc4 ​​values ​​and SA3 structures belonging to C1 or C4, seven drugs were found to be positive for DIIC: amiodarone (CID: 2157), methoxsalen (CID: 4114), dronedarone (CID: 208898), lapatinib (CID: 208908), nitrofurantoin (CID: 6604200), dantrolene (CID: 6914273), and cefuroxime (CID: 5479529). For example, lapatinib (CID: 208908) significantly inhibits ABCB1 and ABCG2 transporters, further impairing lipid excretion, and has been reported to exhibit significant DIIC. Its metabolite, 8-hydroxylapatinib, has also been reported to interact with reactive species such as aldehydes or quinoneimines, and with transporter proteins, resulting in significant toxicity.

[0110] In summary, characterizing drugs using multidimensional features and employing the KNN-GCN method to classify the DIIC clusters to which each drug may belong helps improve the accuracy of DIIC hepatotoxicity prediction. This method not only successfully predicts the DIIC of a compound but also logically infers the structure-pathway-toxicity relationship. These insights can serve as valuable methods for drug development and evaluation. The multidimensional feature model established in this patent, combined with domain representation, provides a new reference anchor pathway analysis model for exploring drug toxicity.

[0111] The above specific examples are merely 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 examples that are varied and modified still fall within the scope of protection 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 using medical subject headings in disease databases. Drugs with known DIIC signatures and their targets were also collected and standardized. (2) Calculation and screening of drug characteristics ① Clustering of DIIC risk genes, and performing gene ontology and pathway enrichment analysis on disease risk genes in the g:Profiler database to obtain disease clusters, and calculating the relative average shortest path length between each disease cluster and the drug target in the PPI network using the network proximity algorithm to characterize the biological network topology of the drug, using Z dc It means that, calculated by formula 1-formula 2, ② SARpy was used to identify the structural alerts of DIIC and calculate the IG value and PR value of each structural alert, where IG represents information gain and is calculated using formulas 3-5; PR represents the positive rate and is calculated using formula 6. The IG value and PR value were further weighted according to formula 7 to obtain a comprehensive value, which was expressed as SA. At the same time, the Mordred descriptor corresponding to Mordred was calculated and preprocessed using the Python package of Mordred. Most drugs were null values, descriptors with variances less than 0.1, and descriptors with Pearson coefficients greater than 0.9 were deleted. Finally, SA and Mordred descriptors jointly characterized the molecular structural characteristics of the drug. H(X)=-p1log2p1-p0log2p0 (Formula 3) IG(T)=H(X)-H(X|T) (Formula 5) SA=0.7PR+10IG (Formula 7) ③ Comprehensively consider the biological network topology and molecular structure characteristics of the drug, and use recursive feature elimination cross-validation 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 built using only the Mordred descriptor. The accuracy, precision, recall, specificity, F1 score, Matthews correlation coefficient, and area under the receiver operating characteristic curve 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, thus forming a domain. Then, drugs within each domain are modeled and evaluated separately. At the same time, a random sample 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, wherein The medical subject headings refer to the standardized terms for intrahepatic cholestasis in the MeSH database; disease databases include OMIM, GeneCards, DisGeNET, and CTD databases; drugs are collected in DILIrank and LiverTox and their tags are checked in metaADEDB, and their targets are collected in DrugBank and DGIdb; Normalization refers to converting Uniprot IDs of all drug targets and disease risk genes into standardized gene names in Uniprot.

3. The method according to claim 1, wherein 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, TNR is true negative rate, Accuracy = (TP + TN) / (TP + TN + FP + FN) (Formula 8) P = Precision = TP / (TP + FP) (Formula 9) Recall = Sf = TP / (TP + FN) (Formula 10) Specificity=SP=TN / (TN+FP) (Formula 11) F1=(2×Precision×Recall) / (Precision+Recall)(Formula 12) 4. The method according to claim 1, wherein 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 layer 1 to layer 1+1 is shown in formula 15.

5. The method according to claim 1, wherein The contribution of each input feature in the model was quantified using the SHAP method. The attribution of each feature of the model was plotted in a SHAP beeswax plot, and the impact of the selected features on the pathogenesis 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.

Images