In this study, only medications with human nephrotoxicity data were included. The nephrotoxic structures were extracted from the Side Effect Resource (SIDER) database (Kuhn et al., 2016). SIDER is a widely used database of adverse drug reactions (ADRs), which contained the information on approved drugs and the ADRs on humans. Herein, we retrieved the entire SIDER database and collected those drugs with nephrotoxicity related ADRs with frequency ≥0.1% in the real world. The corresponding structures of included chemicals were downloaded in smiles format from the PubChem database (Kim et al., 2016). The non-nephrotoxic structures were extracted from Zhang’s work (Zhang et al., 2019). They built the negative drug dataset with the drugs without nephrotoxicity from the SIDER database. The included structures were carefully prepared as follows: (1) removing duplicate substances; (2) keep only the main ingredients in mixtures; and (3) salts were converted to their parent forms.

The sufficiently structural diversity was a key issue for global QSAR models to ensure a reasonable predictive accuracy (Ancuceanu et al., 2019). Principal component analysis (PCA) is a well-known technique for reducing the dimensionality and increasing interpretability, which can solve an eigenvector problem by creating new uncorrelated variables that successively maximize variance (Jolliffe and Cadima, 2016). In this study, the chemical space of data sets was analyzed with the first two principal components of CDK (Chemistry Development Kit) Descriptors. The PCA was performed using SPSS Statistics 26.

The model building was performed on the online chemical database and modeling environment (OCHEM), which is a user friendly web-based platform for automatic and simple QSAR modeling (Sushko et al., 2011). OCHEM supports the typical steps of QSAR modeling, and the models can be published and publicly used on the web (Oprisiu et al., 2013; Cui et al., 2019; Pawar et al., 2019; Cui et al., 2021; Hua et al., 2021; Huang et al., 2021; Ta et al., 2021). Among the many state-of-the-art modeling methods available on OCHEM, we applied five widely used traditional machine learning (ML) approaches and five different deep learning (DL) algorithms. As an application of artificial intelligence (AI), ML has been an effective tool for modeling in computational toxicology. Five highly effective and robust ML approaches were used in this study, including associative neural network (ASNN) (Tetko, 2009), support vector machine (SVM) (Chang and Lin, 2011), C4.5 decision tree (WEKA J48) (Hall et al., 2009), random forest (RFR) (Breiman, 2001), and extreme gradient boosting (XGBoost) (Chen and Guestrin, 2016). DL is an extension of machine learning, and its concept came from the research of artificial neural networks. In this study, we used five different DL approaches, including convolutional neural network fingerprint (CNF) (Tetko et al., 2019), transformer convolutional neural network (TRANSNN) (Karpov et al., 2020), TEXTCNN algorithm available from DeepChem (TEXTCNN) (Wu et al., 2017), Graph Isomorphism Network (GIN) (Capela et al., 2019), and edge attention based multi-relational graph convolutional networks (EAGCNG) (Shang et al., 2018). The detailed descriptions of these ML and DL approaches can be found in the corresponding literature.

The individual parameters for each model algorithm were optimized automatically by the method itself in an inner loop of cross-validation. For instance, the SVM algorithm used libSVM (Chang and Lin, 2011) on OCHEM. There were two important configurable options for this method, including SVM type (ε-SVR and μ-SVR, etc.) and the kernel type (linear, polynomial, radial basis function, sigmoid, etc.). In the OCHEM workflow, classic ε-SVR and radial basis function kernels were used. The other SVM parameters, namely cost C and width of the RBF kernel (γ, g), were optimized using default grid search, and this was performed according to the LibSVM manual (Tetko et al., 2013).

For the machine learning modeling, the molecular descriptors were served as the input of drug structures. We calculated eight different descriptor packages, including Chemaxon descriptors (Chemaxon, 499 descriptors), Fragmentor, GSFrag descriptors (GSFrag, 1,138 descriptors), MORDRED descriptors (MORDRED, 1826 descriptors), PyDescriptor (1,624 descriptors), QNPR descriptors (QNPR), RDKit descriptors (RDKit), and alvaDesc descriptors (5,666 descriptors). The details of these descriptor packages can be learned via http://docs.ochem.eu//display/MAN/Molecular+descriptors.html. The descriptors were filtered with pairwise de-correlation method before the model building. There was no selection bias, since the unsupervised filtering was totally independent.

For the deep learning models, the SMILES string of each compound was served as the input without descriptors.

Consensus modeling is to unify the prediction of unknown samples by multiple individual models to achieve a unified result, and then improve the prediction of the model. By averaging the prediction of individual models, noise can be reduced, thus, and the consensus model can provide better predictive power than most individual models alone (Tapia Garcı’a et al., 2012). Herein, the consensus model was built with simple average of predictions from the best performed individual models.

For QSAR models, it is important to estimate the applicability domain (AD) to determine whether the test compound is suitable for the specific model. In this study, we determined the AD with distance to model (DM) in OCHEM proposed by Sushko (Sushko et al., 2011). The DM assesses the distance from the target compound to the model. The larger DM means the lower applicability for the target compound.

The molecular properties of compounds always can make a significant difference in the toxicity. In this study, the analysis of differences of molecular properties between nephrotoxic and non-nephrotoxic structures was performed, in order to investigate the relevance of these molecular properties with drug induced nephrotoxicity. Several commonly used physical-chemical properties which have been widely adopted in the analysis for other endpoints were calculated and analyzed, including molecular weight (MW), molecular polar surface area (MPSA), AlogP, molecular solubility (LogS), the number of hydrogen bond acceptors (nHBA) and donors (nHBD), the number of rotatable bonds (nRotB), and the number of aromatic rings (nAR). These properties are relative to the molecular size, lipophilicity and solubility, and hydrogen bonding ability and complexity, respectively. The comparison between groups was tested by T-test, and the p value <0.05 was considered to indicate statistical significance. The molecular properties were calculated with PaDEL-Descriptor package (Yap, 2011).

Structural alert (SA), or privileged substructure, was defined as the substructure which can cause the chemicals to become toxic. It has been well accepted in toxicity assessment, because of the direct derivation from mechanistic knowledge. SAs have been commonly used for toxicity assessment of many different endpoints (Claesson and Minidis, 2018; Yang et al., 2018; Wedlake et al., 2020). In this study, we identified the structural alerts for nephrotoxicity by calculating f-score and positive rate of each fragment from Klekota-Roth fingerprint (KRFP, 4,860 bits). The specific substructure should be regarded as a SA if presented more frequently in nephrotoxic drugs than non-nephrotoxic drugs. The positive rate (PR) of a substructure was defined as below: P R = N f r a g m e n t _ p o s i t i v e N f r a g m e n t (1) where N _{fragment_positive} was the number of nephrotoxic drugs containing the fragment, and N _fragment was the total number of drugs containing the fragment.

All the ML and DL models were first internally validated with 5-fold cross validation, and the best performed models were further validated with the external validation set. The structural alerts were assessed with the whole dataset. The classifiers and SAs were evaluated based on the counts of true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN). Several statistical parameters were calculated to evaluate the classifiers, including the total accuracy (Q), sensitivity (SE), specificity (SP), enrichment factor (EF), and Matthews correlation coefficient (MCC), which were calculated with Eqs (2–6). Q = T P + T N T P + F N + T N + F P (2) S E = T P T P + F N (3) S P = T N T N + F P (4) E F = T N / ( T N + F N ) ( T N + F P ) / ( T P + F N + T N + F P ) (5) M C C = T P ∗ T N − F P ∗ F N ( T P + F P ) ( T P + F N ) ( T N + F P ) ( T N + F N ) (6)

Additionally, the receiver operating characteristic (ROC) curve was also plotted for the QSAR models, and the values of area under the ROC curve (AUC) were provided, too.

In the study, 565 diverse structures were kept after preparation, including 287 nephrotoxic drugs and 278 non-nephrotoxic drugs. As shown in Table 1, the whole data set was randomly divided into a training set with 456 chemicals (232 nephrotoxic and 224 non-nephrotoxic) and an external validation set with 109 chemicals (55 nephrotoxic and 54 non-nephrotoxic). The structures of the drugs can be found in Supplementary Table S1.

The diversity of structures is a crucial factor for the applicability of global models. Thus, we performed the principal component analysis (PCA) based on CDK descriptors to analyze the chemical space of the included compounds. PCA can simplify the complexity in high-dimensional data while retaining trends and patterns by transforming the data into fewer dimensions, which act as summaries of features (Ringnér, 2008). The first two principal components were kept and used for the definition of the chemical space of compound in the training and validation sets in this study. As shown in Figure 1, the distribution scatter diagram illustrated that the data sets shared a similar chemical space. In addition, the Tanimoto similarity index (TSI) was also calculated based on the ECFP-4 fingerprint to evaluate similarities among the structures. The average value of the entire data set was 0.13, which indicated an evidently chemical diversity of the entire data set.

Combined with five different ML algorithms and eight molecule descriptor packages, 40 ML models were developed. Meanwhile, five DL models were developed with different algorithms using chemical SMILES-string as input. The performances of ML and DL models on 5-fold cross validation were shown in Table 2. The prediction accuracy (Q) of the models ranged from 58.54 to 73.90%. Among them, three models performed much better on 5-fold cross validation than others, including a DL model (CNF) and two ML models (XGBoost_QNPR and RFR_QNPR). These models all provide good predictive ability with Q value ≥ 70.00% and AUC value ≥0.80.

The consensus approach has been proved able to improve the accuracy of models. In this study, a consensus model was developed based on the three best models. The consensus model performed much better than the individual models. As shown in Table 2, it provided the Q value of 75.88% and the AUC value 0.83; the values of SP, SE, EF, and MCC were 72.84, 79.01, 1.50, and 0.52%, respectively.

Due to the complete independence from the model training, the external validation set could be well used for evaluating the predictive ability of models objectively. As shown in Table 3; Figure 2, the best performed individual models also achieved good predictive results on external validation. The model developed with RFR algorithm and QNPR descriptors performed best with Q value 87.16% and AUC value 0.91. The DL model developed with CNF algorithm also provided good predictive ability, with Q value 83.49%, and AUC value 0.89. The consensus model did not perform better than RFR_QNPR model on most of the statistical parameters, except for AUC value (0.93). It provided a Q value of 86.24% and MCC value of 0.82 on external validation, and the values of SE, SP, and EF were 85.45, 87.04, and 1.72%, respectively.

The results of ML model building and validation suggested that ML models developed with the QNPR descriptor performed much better than others. The QNPR descriptors were derived directly from the chemical SMILES strings. For each structure either canonical SMILES or IUPAC name would be split into fragments of a specified length determined by the configuration (Thormann et al., 2007). For the different ML algorithms, XGBoost and RFR performed better than others. XGBoost is a scalable end-to-end tree boosting system, which has been widely used in many application scenarios. It supports multi-thread computation and uses regularization enhancement technology to reduce over-fitting, so as to ensure the robustness of the model. Meanwhile, it has the advantages of flexibility, fast calculation speed, and good robustness. Therefore, XGBoost is not easy to be disturbed by outliers. It can achieve state-of-the-art results on many machine learning challenges (Chen and Guestrin, 2016). RFR is an ensemble learning method by constructing a multitude of decision trees (DT) at training time (Breiman, 2001). It can help to improve the accuracy by reducing overfitting in decision trees and automating missing values present in the data. CNF contributed to the best DL model. CNF is one of the state-of-the-art deep learning methods, which is based on ideas of text processing. It could achieve the high prediction accuracy due to the augmentation technique employed during both training and inference steps.

The DL model in our study did not show significant better predictive power compared with the ML model. However, considering that the DL models did not require additional molecular description, as long as SMILES are provided, the DL algorithms still showed obvious advantages over ML algorithms.

The molecular physical-chemical properties could provide useful information for biological activities. Herein, we analyzed the differences of several commonly used physical-chemical properties between nephrotoxic and non-nephrotoxic drugs. The distributions of these descriptors for nephrotoxic and non-nephrotoxic drugs can be seen in Figure 3.

The characteristics MW and MPSA can simply assess the size and complexity of compounds. For the entire data set in the study, the values of MW were distributed between 59.04 and 4,491.88, with a mean of 416.86. The mean value was 459.71 for nephrotoxic drugs and 372.63 for non-nephrotoxic drugs. The difference between the mean MW of nephrotoxic and non-nephrotoxic drugs was not significantly different (p < 0.01). The values of MPSA were distributed from 0 to 1902.88, with a mean of 119.80 for the entire dataset. The mean value was 137.52 for nephrotoxic drugs and 101.51 for non-nephrotoxic drugs (p < 0.01). These results suggested that it may be significantly different in structure size and polar surface area between nephrotoxic and non-nephrotoxic compounds.

The AlogP is commonly used to represent the lipophilicity of compounds. For the entire data set, the values of AlogP ranged from −32.81 to 10.25, with a mean of 1.48. The mean value of AlogP was 1.34 for nephrotoxic chemicals and 1.64 for non-nephrotoxic drugs. As can be seen in Figure 3, the distributions were not significantly different between nephrotoxic and non-nephrotoxic drugs with p-value 0.34. The result indicated that chemical lipophilicity may be weakly correlated with drug induced nephrotoxicity.

Chemical hydrogen bonding ability also is an important character for its activity and toxicity, and it was usually in terms of nHBA and nHBD. In this data set, mean values of nHBA for nephrotoxic and non-nephrotoxic drugs were 7.06 and 5.41, respectively, and the difference was significant (p < 0.01). Meanwhile, the nephrotoxic and non-nephrotoxic drugs had mean values of nHBD with 3.29 and 2.47, respectively, and the difference is not significant (p = 0.07). The data indicated that nHBA was obviously associated with drug induced nephrotoxicity while nHBD was not.

LogS is an estimation of molecular solubility in water. For the entire data set, the values of LogS were distributed from −35.86 to 2.45, and the mean value was −4.52. For nephrotoxic structures, the mean value of LogS was −4.81, and it is −4.21 for non-nephrotoxic drugs. The difference also proved statistical significance (p = 0.03). This result demonstrated that there was difference of molecular solubility between nephrotoxic and non-nephrotoxic drugs.

As shown in Figure 3, drug induced nephrotoxicity was also obviously associated with nRotB (the mean values were 8.05 for nephrotoxic chemicals, and 5.78 for non-nephrotoxic drugs, and with p-value 0.01) and nAR (the mean values were 1.70 for nephrotoxic chemicals and 1.08 for non-nephrotoxic drugs, with significant difference (p < 0.01).

Through the analysis of molecular properties, several physical-chemical properties have obvious differentiating effect on drug induced nephrotoxicity. Nevertheless, nephrotoxicity is a complex endpoint yet. It’s not easy to explain the mechanism of drug induced nephrotoxicity with individual simple chemical descriptors.

The structural alerts (SA) responsible for nephrotoxicity were identified using f-score and positive rate analysis of each fragment from KRFP fingerprint. Only the fragments existing in 6 or more drugs were kept. The fragments with f-score ≥0.005 and positive rate ≥0.75 were identified. Finally, 87 representative fragments were filtered and listed in Supplementary Table S2. Among them, 16 substructures presented in nephrotoxic active chemicals only, which covered 76 nephrotoxic drugs. Details of each fragment and the representative structures were shown in Table 4.

We analyzed the structural characteristics of drugs in the entire data set, and 75.04% structures (424/565) were correctly classified. Among the 290 drug structures contained at least one identified substructure, 218 structures were true nephrotoxic. The classification accuracy was 75.17%. For the drug structures that did not contain any identified substructure, 206 of 275 drugs (74.91%) were true non-nephrotoxic. The frequencies of these substructures were much higher in nephrotoxic drugs than non-nephrotoxic drugs, and showed a good ability to distinguish nephrotoxic drugs in the whole data set. To a certain degree, these fragments could be considered as the structural alerts responsible for nephrotoxicity. If a structure contains one or more structural alerts, it is more likely to be nephrotoxic than non-nephrotoxic.

Fluorine is widely used in medicinal chemistry to improve a molecule’s potency and permeability. However, the metabolism of fluorinated compounds may produce fluoride and other toxic metabolites (Pan, 2019). As shown in Table 4, several selected structural alerts contained phenyl fluoride (Nos 3, 8, 9, 11, 12, and 14). Since the kidney is a main target organ of mammalian fluoride systemic exposure and renal toxicity can occur after acute and chronic fluoride intoxication (Quadri et al., 2016). The kidney plays a vital role in fluoride metabolism, since 50–80% of fluoride is removed via urinary excretion. Fluorinated compounds can be toxic to the kidneys in a number of ways. It has widely been reported that fluorinated compounds can increase the generation of reactive oxygen species (ROS) and free radicals, cause extensive oxidative stress and excessive lipid peroxidation, and reduce antioxidant enzyme activities in vivo or in vitro (Quadri et al., 2016). Thus, numerous renal structural, ultrastructural, and functional may be changed after receiving increased amounts of fluorinated compounds exposure. Certainly, this is not meant to sound the alarm on all fluorinated compounds, we suggest raising the awareness of common drug instability and metabolism issues leading to defluorination, as well as the resulting reactive/toxic metabolite (Pan, 2019). Polyamines and their derivatives (Nos 2 and 15) widely existed in aminoglycosides and the other nephrotoxic drugs. Polyamines catabolism can be stimulated through oxidation, which will lead to the generation of ROS and low abundance of free polyamines with antioxidant capacity (Pegg, 2013; Murray Stewart et al., 2018). Consequently, this process is associated with a number of pathologies, including cancer, neurological disorders, and kidney dysfunction. Recent studies have also suggested a role for polyamines derivatives in the p53-mediated ferroptotic response to ROS stress. It is an iron-dependent and nonapoptotic mode of cell death and has been associated with drug-induced nephrotoxicity recently (Ou et al., 2016; Murray Stewart et al., 2018). Benzimidazoles derivatives were always proposed as new bioreductive prodrugs with the potential anticancer activity, due to their effect on the DNA destruction and growth inhibition into selected tumor cell lines (Błaszczak-Świątkiewicz et al., 2014). Besides, benzimidazoles derivatives also can induce concurrent apoptotic and pyroptotic cell death (Ren et al., 2021). However, when these processes take place in kidney cells, kidney damage can occur. Toluene-contained structures (Nos 1 and 7) were also selected as structural alerts for drug induced rhabdomyolysis (Cui et al., 2019), which always associated to oliguric renal failure. In the setting of toluene intoxication, electrolyte disturbances may play important roles on causing rhabdomyolysis (Camara-Lemarroy et al., 2015). Toluene has also been associated with direct induction of acute tubular necrosis and acute oliguric renal failure. Among the SAs we identified, there were several fragments that also could cause the formation of ROS (Nos 6 and 16, etc.), which may contribute to the nephrotoxicity.

In the present study, the SAs were identified with the f-score and frequency analysis of defined substructures. These methods have been widely used for the SA discovery for many other endpoints. In fact, these methods have some shortcomings. They are not able to characterize the spatial arrangement of identified substructures, and they cannot make a good distinction when more than one SA presented in the same structure. In spite of this, these structural alerts were able to well distinguish chemical structures with renal toxicity, and they can help to understand the specific fragments which lead to nephrotoxicity. Therefore, these SAs should be severed as a useful tool to visually evaluate the nephrotoxicity of chemicals.

For ease of use, the QSAR models were made available at OCHEM. The consensus model could be accessed via https://ochem.eu/article/140251. The three best individual models (RFR_QNPR, XGBOOST_QNPR, and CNF) were also available with the corresponding model IDs. Users can predict chemical nephrotoxicity by using the “Apply the model to new compounds” link. In addition, the data sets for modeling could be downloaded by using the “Export this basket” link.

The structural alerts responsible for nephrotoxicity have been integrated as part of our Web server SApredictor, which is a structural alert based expert system for drug toxicity prediction and freely available at http://www.sapredictor.cn. With the help of SApredictor, people can quickly evaluate whether the query chemicals are nephrotoxic, and the specific structural fragments that lead to the nephrotoxicity of the compounds will be intuitively shown to provide valuable reference for the modification of the structures.