{"paper_id":"b65da321-cae8-4a85-9db8-e5df40293430","body_text":"Gonadotropin-releasing hormone (GnRH) helps regulate the hypothalamic–pituitary–gonadal axis, which governs reproduction and sexual development. GnRH, produced by neurosecretory cells in the hypothalamus, binds to the GnRH receptor (GnRHR) expressed on gonadotroph cells in the anterior pituitary gland and acts in a pulsatile manner [ 1 ]. Consequently, the anterior pituitary gland secretes gonadotropins such as follicle-stimulating hormone and luteinizing hormone, which both act on the ovaries or testes and promote estrogen and testosterone production [ 2 ]. GnRH secretion is suppressed during childhood, but it increases at puberty to induce sexual development through increased gonadotropin and gonadal steroid production [ 3 ]. GnRHR signaling dysregulation can cause clinically important reproductive disorders, including polycystic ovary syndrome [ 4 ] and ovulatory dysfunction [ 5 ].\nGnRHR agonists are therapeutic agents that are widely used for inducing ovulation in infertility treatment [ 6 ] and for treating prostate cancer [ 7 ], breast cancer [ 8 ], and endometriosis [ 9 ]. However, these drugs have been associated with adverse effects, including amenorrhea, hot flashes, and reduced bone mineral density [ 10 ]. More recently, they have been reported to be potentially associated with more serious adverse events, including interstitial lung disease [ 11 ] and liver disease [ 12 ]. Furthermore, a large-scale analysis using the Food and Drug Administration (FDA) Adverse-Event Reporting System (FAERS) suggested associations with reproductive system and breast disorders [ 13 ]. Nevertheless, previous studies mainly focused on specific drugs, and the overall adverse-event profile associated with GnRHR agonist activity remains insufficiently characterized.\nThe Toxicology in the 21st Century (Tox21) program was launched through a collaboration between the U.S. Environmental Protection Agency, the U.S. National Institutes of Health, and the U.S. FDA to promote the transition from conventional animal-dependent toxicity testing to in vitro assessment [ 14 ]. As part of this program, Yang et al. conducted comprehensive quantitative high-throughput screening of the Tox21 compound library against human and mouse GnRHRs in 2024; compounds with GnRHR agonist activity were identified from a broad chemical space comprising not only approved drugs but also environmental chemicals [ 15 ].\nLarge-scale in vitro datasets such as Tox21 provide an important foundation for constructing quantitative structure–activity relationship (QSAR) models. QSAR models predict biological activity from molecular structure, with recent studies integrating QSAR-based toxicity prediction with FAERS data [ 16 ]. However, the validity of the evaluation design must be given careful attention when assessing the QSAR model performance. Calle et al. reported that nested cross-validation reduces the performance estimate variance compared with single-split evaluation and avoids test-information leakage by tuning hyperparameters within an inner loop [ 17 ].\nIn this study, we aimed to clarify the potential impact of GnRHR agonist activity on drug safety. We developed a prediction model based on compound data for GnRHR agonist activity provided by the Tox21 program and evaluated its performance using nested cross-validation. The model was applied to drugs listed in the FAERS database, and trends in adverse events associated with predicted GnRHR agonist activity were then analyzed.\n\nThe four models, namely, RandomForest, LightGBM, XGBoost, and BalancedRandomForest, were compared using the mean values of evaluation metrics obtained in the outer loop of nested cross-validation. BalancedRandomForest achieved the best performance in terms of ROC-AUC, ACC, MCC, and SP, with an ROC-AUC of 0.810 ( Table 1 ), and was therefore selected for subsequent drug screening. Confidence intervals for ROC-AUC and PR-AUC based on out-of-fold predictions from the five outer folds are provided in  Supplementary Table S1 . A sensitivity analysis excluding PubChem inconclusive compounds yielded improved predictive performance across all algorithms; however, the ranking of models based on ROC-AUC remained unchanged, with BalancedRandomForest maintaining the highest performance ( Supplementary Table S2 ). In addition, the final model retrained using all data was employed in the analysis. Threshold-performance analysis across multiple uncertainty intervals demonstrated that the 0.3–0.7 range provided a favorable balance between applicability-domain coverage and prediction reliability ( Supplementary Table S3 ). Therefore, compounds within this range were excluded from the applicability domain, while those with scores ≥0.7 were regarded as predicted positives.\nFeature importance was calculated by SHAP analysis using the final model constructed with BalancedRandomForest and all data used for model construction. Descriptors such as fr_Al_OH, SsOH, MAXsOH, PNSA3, VSA_EState3, BCUT2D_LOGLPLOW, GATS5i, GATS4pe, and MATS4Z, which represent solvent interaction properties; MinPartialCharge, MaxAbsPartialCharge, ATSC7c, and MATS7c, which represent local electronic properties; and Xc-3dv, Xc-4dv, Xp-2dv, Xpc-5dv, Chi0v, VE3_Dzse, and ETA_beta_s, which represent skeletal complexity, were the leading features ( Figure 1 ). Higher values of fr_Al_OH, SsOH, MAXsOH, and all skeletal complexity descriptors and lower values of PNSA3, BCUT2D_LOGLPLOW, GATS4pe, ATSC7c, and MATS7c increased the likelihood that a compound would be predicted as a GnRHR agonist.\nIn total, 5523 FAERS-listed drugs for which SMILES strings could be obtained served as candidates for GnRHR agonist prediction. When the applicability domain was defined using the minimum and maximum scores of PC1 and PC2 from PCA based on the Tox21 data, 5177 of these drugs were considered within the applicability domain ( Figure 2 ). When drugs with an intermediate normalized score for the active class (0.3–0.7) from the final BalancedRandomForest model were considered outside the applicability domain, 1191 drugs remained within the applicability domain, of which 367 were predicted to be positive for GnRHR agonist activity ( Supplementary Table S4 ). As a supplementary validation of the applicability domain, a distance-based k-nearest neighbors (kNN) approach was also evaluated. Using the kNN criterion, 4093 FAERS-listed drugs were classified as being within the applicability domain. Among these compounds, 157 were predicted to have GnRHR agonist activity ( Supplementary Table S4 ). Although the kNN-based approach yielded a more restrictive applicability domain than the PCA- and score-based approach, a substantial number of compounds predicted to have GnRHR agonist activity remained within the applicability domain, supporting the robustness of the overall findings.\nTo characterize the 367 drugs predicted to have GnRHR agonist activity, we classified the drugs according to the first level of the Anatomical Therapeutic Chemical (ATC) classification system. Among FAERS-listed drugs, 8066 drugs were assigned ATC codes, and 193 drugs were both predicted to exhibit GnRHR agonist activity and assigned ATC codes. We compared the distribution of predicted GnRHR agonists by ATC first-level category with that of all ATC-classified, FAERS-listed drugs, and the proportion in each ATC category among these FAERS-listed drugs was set to 100% ( Figure 3 ). Drugs predicted to have GnRHR agonist activity were overrepresented among systemic hormonal preparations, excluding sex hormones and insulins (H), whereas drugs acting on the musculoskeletal system (M) and nervous system (N) were underrepresented.\nFirst, the FAERS DRUG table (127,228,343 rows) and REAC table (54,645,478 rows) were merged to construct a case-based database. In this database, the number of cases after duplicate removal was used as the analytical unit, and for each adverse event, 2 × 2 contingency tables were generated according to the presence/absence of exposure to drugs predicted to be positive for GnRHR agonist activity and the presence/absence of the adverse event. From 18,328,780 cases, cases involving the 367 drugs predicted to be GnRHR agonists were extracted. A total of 330 unique PTs, corresponding to 531 PT-SOC assignments, met the following criteria: at least 1000 reports, ln(ROR) > 1, and  q  < 0.05 ( Supplementary Table S5 ). When adverse events with at least 1000 reports were summarized by MedDRA System Organ Class (SOC) ( Table 2 ), significant adverse events occurred most frequently in respiratory, thoracic and mediastinal disorders ( Supplementary Figure S1 ), infections and infestations ( Supplementary Figure S2 ), and hepatobiliary disorders ( Supplementary Figure S3 ).\n\nUsing GnRHR agonist activity data from the Tox21 database, we developed a machine learning model based on 1761 molecular descriptors. In nested cross-validation, the model achieved a mean ROC-AUC of 0.810, indicating a moderate ability to rank active compounds above inactive compounds. However, the relatively low PR-AUC and MCC indicated that the model’s ability to make definitive positive predictions for individual compounds was limited, particularly given the substantial class imbalance. Accordingly, the model was used to prioritize candidate drugs and generate hypotheses through subsequent FAERS adverse event analysis, rather than to confirm the GnRHR agonist activity of individual drugs. The prediction model’s applicability domain was evaluated from two perspectives: (i) distribution in the chemical space based on PCA and (ii) distribution of normalized model scores. In the PCA-based assessment, the chemical applicability domain was defined according to the range of Tox21 compounds in the PC1–PC2 plane. Many compounds predicted to have high GnRHR agonist activity were distributed within or near this region. This result supports the model’s external validity because predictions were made consistently within the chemical space represented by the training data [ 18 ]. However, the variance explained by PC1 and PC2 in this dataset was relatively low (23.5% and 6.7%, respectively); thus, PCA alone may be insufficient to strictly define the applicability domain. Therefore, the distribution of normalized scores was used as a complementary criterion, and compounds with normalized active-class scores within the intermediate range of 0.3–0.7 were considered outside the applicability domain because of high model uncertainty. In classification models, score-based applicability domain definitions reportedly improve predictive performance more effectively than definitions based solely on distance information, including PCA [ 19 ]. Thus, the criteria used in the present study align with previous findings. By defining the applicability domain using both chemical space and normalized model scores, the present model enabled us to select more reliable candidate compounds while reducing excessive extrapolation.\nIn the prediction model, higher values of descriptors related to the number of hydroxyl groups, including fr_Al_OH, SsOH, and MAXsOH, as well as descriptors reflecting skeletal complexity, increased the likelihood of being predicted as a GnRHR agonist. Conversely, lower values of PNSA3 (representing a polar and negatively contributing molecular surface area), BCUT2D_LOGLPLOW (reflecting lipophilicity), GATS4pe (reflecting pi-electronic properties), ATSC7c (reflecting the overall spread of charge), and MATS7c (reflecting charge regularity) were associated with a higher probability of being predicted as a GnRHR agonist. GnRHR contains many polar amino acids and depends on hydrogen bonding and electrostatic interactions [ 20 ]. It is also regarded as a distinctive G protein-coupled receptor that selectively recognizes a U-shaped folded GnRH structure in the ligand-binding domain instead of a disordered peptide conformation [ 21 , 22 ]. Furthermore, GnRH acts as a hydrophilic peptide enriched with polar amino acids, with many GnRH analogs containing multiple hydroxyl groups and forming sterically constrained structures through D-amino acid introduction or cyclization. These observations are broadly compatible with the molecular characteristics of GnRHR agonists identified in our study and suggest that the prediction model captures, at least in part, the molecular properties of GnRHR agonists.\nDrugs predicted to be GnRHR agonists were enriched among systemic hormonal preparations excluding sex hormones and insulins (H), consistent with the representative GnRHR agonists gonadorelin and nafarelin classified as GnRHs (H01CA). Similarly, antineoplastic and immunomodulating agents (L)—the next most frequent category—include buserelin, leuprorelin, goserelin, triptorelin, and histrelin, which are classified as GnRH analogs (L02AE). In vitro GnRHR agonist activity has also been confirmed in the Tox21 database for adrenergic drugs (e.g., racepinephrine and epinephrine), cholinergic drugs (e.g., acetylcholine), and bile acids (e.g., chenodeoxycholic acid). These compounds were predicted by the present model to have GnRHR activity, given that their normalized scores were above 0.95. Therefore, through this analytical framework, GnRHR-associated adverse-event risks may be evaluable not only for established GnRHR agonists but also for drugs with potential GnRHR agonist activity.\nDrugs predicted to be GnRHR agonists were frequently associated with respiratory adverse events in this study. This finding aligns with previous reports and existing evidence linking androgen deficiency to effects on lung tissue.\nSeveral cohort and database studies have reported that GnRHR agonist therapy is associated with increased risks of interstitial lung disease [ 11 ], community-acquired pneumonia [ 23 ], and radiation pneumonitis [ 24 ].\nOjeda et al. assessed morphological and biochemical changes in the lungs of male rats receiving surgical or pharmacological castration. They found that androgen deficiency induced alveolar septal destruction, emphysema-like changes, and alterations in the phospholipid composition of pulmonary surfactant [ 25 ]. Rats treated with antiandrogen flutamide also demonstrated similar morphological lung changes. Therefore, androgens are potentially essential for maintaining normal lung morphology and pulmonary surfactant metabolism.\nTaken together, androgen deficiency induced by GnRHR agonists may compromise the structural integrity of lung tissue through morphological and biochemical changes in pulmonary cell membranes. This may increase susceptibility to respiratory diseases, including pneumonia.\nThe high number of reported infection-related adverse events aligns with a series of evidence linking androgens to immune function.\nAndrogen receptors are expressed in many immune cells, including neutrophils [ 26 ], and androgens promote granulocyte progenitor proliferation and increase neutrophil production [ 27 ]. Using androgen receptor–knockout mice, Chuang et al. showed that androgen receptor deficiency was associated with neutropenia and increased vulnerability to microbial infection. They also showed that exogenous androgen administration restored neutrophil counts in castrated mice [ 28 ].\nAccording to Ojeda et al., androgen deficiency induces changes in the phospholipid composition of pulmonary surfactant and morphological changes in lung tissue [ 25 ], potentially impairing airway defense mechanisms. Hicks et al. also showed that androgen deprivation therapy (ADT) significantly increased the risk of hospitalization for community-acquired pneumonia [ 23 ]. Thus, the synergistic effects of immune dysfunction and pulmonary tissue changes may increase infection risk among GnRHR agonist users.\nHepatobiliary adverse events were also frequently reported for drugs predicted to be GnRHR agonists in this study, consistent with the known metabolic adverse effects of ADT.\nIn a large population-based study, Gild et al. showed that ADT, including GnRHR agonists, was associated with liver disease diagnoses such as cirrhosis, hepatic necrosis, and nonalcoholic fatty liver disease in men with prostate cancer [ 12 ]. Saylor et al. also reported that GnRHR agonist therapy significantly increased biliary disease incidence and that the risk increased with longer duration of use, suggesting a dose–response relationship [ 29 ]. In addition, bilateral orchiectomy was not associated with a significant increase in biliary disease risk, suggesting that mechanisms specific to GnRHR agonists contribute to biliary disease.\nMoreover, GnRHR agonists induce metabolic abnormalities, including increased fat mass, particularly abdominal fat accumulation, elevated serum triglycerides, and increased insulin resistance [ 30 ]. According to a metabolomic analysis, most bile acids and their metabolites increase within 3 months after the initiation of GnRHR agonist therapy [ 31 ]. Additionally, hormonal changes can induce gallbladder dysmotility [ 32 ]. These factors may have contributed to the increased risk of developing hepatobiliary diseases observed in the present study.\nIn the present study, the proportion of reported adverse events classified as reproductive system and breast disorders was relatively low among cases with GnRHR-active drug exposure ( Supplementary Figure S4 ). However, this result may reflect changes in the composition of adverse-event reports related to underlying indications rather than the drugs’ safety profile. Given that GnRHR agonists are mainly used to treat breast cancer [ 8 ] and prostate cancer [ 7 ], reports related to these diseases may have influenced the overall reporting structure, making the proportion of reproductive system and breast disorders appear relatively low.\nConversely, ovarian hyperstimulation syndrome (ROR = 13.48,  q  < 3.41 × 10 −322 ) and polycystic ovaries (ROR = 3.17,  q  = 4.06 × 10 −78 ), which are known adverse events linked to GnRHR agonists, were detected as significant signals. These findings align with previous reports on ovarian stimulation and effects on reproductive endocrine function induced by GnRHR agonists [ 1 , 2 , 3 , 6 ], supporting the validity of the present analysis.\nTherefore, the lower proportion of reported reproductive system and breast disorders may reflect the influence of underlying indications and changes in the reporting structure rather than a reduced risk of reproductive system-related adverse events.\nThis study has several limitations that need to be considered when interpreting and generalizing the findings. First, machine learning models based on QSAR may be biased by the training data’s structural distribution, and their predictive performance for structurally distinct unknown compounds may be limited. Furthermore, some molecular descriptors are expected to be highly correlated, potentially affecting the model’s interpretability and reproducibility [ 33 ]. In addition, the 3D molecular descriptors were calculated from MMFF-minimized conformers, which may not represent the actual receptor-bound conformations of ligands.\nSecond, given that Tox21 activity data are based on in vitro assays, the direct relationship of these data with in vivo pharmacological activity or actual adverse effects remains unclear. In addition, the inclusion of PubChem inconclusive compounds in the negative class may have introduced label ambiguity and reduced class separability, which could have affected the model’s generalizability and classification performance. These limitations affect how the prediction results are interpreted clinically. Third, the model was applied to 1191 drugs for which SMILES strings were available and were within the applicability domain of the Tox21-based prediction model. Consequently, other drugs were excluded from the prediction results, and the influence of these excluded compounds on the overall evaluation remains unknown.\nFourth, the present analysis was conducted at the MedDRA SOC level and may not fully capture variation among individual adverse events. Pharmacovigilance research may exhibit a masking effect, wherein a large number of other adverse-event reports obscure specific adverse-event signals [ 34 , 35 ]. Indeed, in our analysis, the proportion of reported reproductive system and breast disorders was low at the SOC level, while significant signals were detected for individual PTs, including ovarian hyperstimulation syndrome and polycystic ovaries. Therefore, SOC-level aggregation may not fully capture detailed changes in adverse-event profiles.\nFifth, spontaneous reporting databases record only cases with reported adverse events and without information on all patients exposed to a drug [ 36 ]. Thus, the incidence of adverse events cannot be calculated directly, and indicators such as ROR must be used for risk evaluation, thereby potentially introducing bias. Sixth, reporting bias must be considered. Some adverse events, particularly mild ones, or newly emerging risks may be underreported, whereas serious adverse events or known risks may be overreported or underreported, likely affecting the results [ 36 , 37 ]. Seventh, indication bias and channeling bias should be considered when interpreting the results. The observed disproportionality signals may be influenced by the underlying diseases, patient characteristics, and clinical contexts for which the drugs were prescribed [ 34 ]. Furthermore, patients may be systematically selected for GnRHR agonist treatment on the basis of clinical characteristics, including disease severity, comorbidities, concomitant therapies, and monitoring intensity, potentially affecting the observed reporting patterns.\nFinally, when multiple drugs are used concomitantly, identifying which drug caused an adverse event is difficult [ 36 , 38 ]. Our analysis included drugs recorded as PS, secondary suspect, concomitant, or interacting, thereby potentially leading to increased FP findings. Therefore, causal relationships between drugs and adverse events cannot be established from this study alone. Further research that accounts for these factors is warranted.\n\nWe collected data on GnRHR agonists from the Tox21 database (AID: 1920067;  https://pubchem.ncbi.nlm.nih.gov/bioassay/1920067 ; accessed on 6 January 2025) to comprehensively analyze compounds with GnRHR agonist activity. Based on high-throughput screening, Tox21 is a large-scale toxicity database developed for an efficient and predictive assessment of chemical toxicity [ 14 ]. From this database, compounds with Simplified Molecular Input Line Entry System (SMILES) strings and activity scores were included, whereas those with duplicate SMILES strings and without carbon atoms were excluded. Subsequently, 7133 compounds were included in the analysis. For compounds with duplicate SMILES strings, the median PubChemScore was calculated. According to PubChem criteria, compounds with PubChemScore ≥40 were considered positive (label: 1), whereas those with PubChemScore <40 were considered negative (label: 0). The negative class included 252 compounds designated as inconclusive by PubChem (activity score ≥1 and <40). Ultimately, the final dataset comprised 7071 negative and 62 positive compounds.\nRDKit, an open-source cheminformatics toolkit, was used for cleaning and standardizing SMILES strings [ 39 ] by removing salts, counterions, and fragments and adjusting protonation states by neutralization. Considering the computational burden associated with a multitude of candidate compounds, optimal three-dimensional (3D) structures were generated using an efficiency-oriented heuristic procedure. First, chemical structures were generated from the SMILES strings, and explicit hydrogen atoms were added. Next, up to 200 3D conformers were randomly generated. Each conformer was energy-minimized using the Merck molecular force field (MMFF), and the one with the lowest energy was selected. If this process exceeded 60 s, conformers were generated using the Experimental-Torsion Knowledge Distance Geometry method [ 40 ], followed by energy minimization with MMFF [ 41 ]. The selected conformer was then converted into a structure-data file (SDF) format.\nFor each compound, molecular descriptors were calculated using the Python libraries RDKit and mordred-community, an actively maintained fork of Mordred [ 42 , 43 ]. Two-dimensional and 3D descriptors, including RDKit descriptors (221 dimensions) and Mordred descriptors (1826 dimensions), were obtained. From a total of 2047 descriptors, perfectly correlated descriptors and those with zero variance were removed, leaving 1761 descriptors for model construction. These preprocessed descriptors were used for machine learning model development without further scaling. For PCA, the descriptor values were standardized to zero mean and unit variance prior to analysis.\nClassification models were developed to predict GnRHR agonist activity. Four algorithms were compared: RandomForest [ 44 ], LightGBM [ 45 ], XGBoost [ 46 ], and BalancedRandomForest [ 47 ]. The 1761 molecular descriptors calculated using RDKit and Mordred served as the explanatory variables, whereas the presence/absence of GnRHR agonist activity, as defined by PubChemScore, was the objective variable.\nAs an external-validation-like framework, model performance was evaluated using nested cross-validation to estimate generalization performance while reducing overfitting [ 17 ]. Specifically, stratified  k -fold cross-validation was performed using five outer and four inner folds. We used the outer folds for evaluating generalization performance and the inner folds for optimizing hyperparameters. Hyperparameters were optimized using Bayesian optimization with Optuna [ 48 ]. The optimization process included class imbalance handling parameters, where the sample weight assigned to positive samples was varied from 1 to 10 times that assigned to negative samples for all algorithms. For BalancedRandomForest, the sampling_strategy parameter was additionally optimized. For each algorithm, 100 Optuna trials were conducted, and the parameter set with the highest area under the receiver operating characteristic curve (ROC-AUC) was selected. In comparing models, the algorithm with the highest mean ROC-AUC across the five outer folds was selected for the final model. The final model was retrained using all data. A total of 100 hyperparameter patterns were explored with Optuna; the parameter set that achieved the highest ROC-AUC in fivefold cross-validation was selected as the final setting. Feature importance was then calculated by SHapley Additive exPlanations (SHAP) analysis using the final model and all data used for model construction [ 49 ].\nThe classification models’ predictive performance was evaluated according to the number of true positives (TP; positive compounds correctly predicted as positive), true negatives (TN; negative compounds correctly predicted as negative), false negatives (FN; positive compounds misclassified as negative), and false positives (FP; negative compounds misclassified as positive). The following seven metrics were used for the evaluation.\nSensitivity (SE): SE is the proportion of actual positive samples correctly predicted as positive. (1) S E = T P T P + F N\nSpecificity (SP): SP is the proportion of actual negative samples correctly predicted as negative. (2) S P = T N T N + F P\nAccuracy (ACC): ACC is the number of correctly predicted samples divided by the total number of samples. (3) A C C = T P + T N T P + T N + F P + F N\nBalanced accuracy (BA): BA is the mean of SE and SP. (4) B A = 1 2 S E + S P\nMatthews correlation coefficient (MCC): MCC is a metric for evaluating classification performance in imbalanced datasets [ 50 ]. (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 )\nROC-AUC: ROC-AUC is calculated as the area under the ROC curve based on SE and 1 − SP across thresholds [ 51 ].\nAverage precision (AP): AP is calculated as the area under the precision–recall curve based on precision [TP/(TP + FP)] and recall (SE) across thresholds [ 52 ].\nThe optimal cutoff value for TP, FN, TN, and FP was determined using Youden’s index, which maximizes SE + SP − 1 [ 53 ]. Given that the final prediction model had a model-specific cutoff, the model output probabilities were standardized for drug prediction using the following equation: (6) X n = X u − log c 2 \nwhere  X n  represents the normalized score,  X u  represents the unnormalized model output probability, and  c  represents the cutoff value determined for the final prediction model. The exponent term (−log c 2) was chosen so that the model-specific cutoff value was mapped to a normalized score of 0.5. Consequently, prediction scores from models with different cutoff values could be interpreted on a common scale, with values above 0.5 indicating predictions above the model decision threshold and values below 0.5 indicating predictions below the threshold. This normalization facilitated consistent interpretation and comparison of prediction scores for FAERS-listed drugs. It is intended for score rescaling rather than probability calibration.\nThe FAERS database was accessed from the official FDA website [ 54 ], and data recorded from the first quarter of 2004 through the third quarter of 2024 were used. Adverse events in the reaction (REAC) table are registered according to MedDRA Preferred Terms (PTs), and drug names in the DRUG table were curated by Intage Healthcare Co., Ltd. (Tokyo, Japan). Each case is assigned a PRIMARYID, allowing the integration of data tables. In the analysis, the DRUG and REAC tables were merged by PRIMARYID. When the same drug and the same adverse event were recorded multiple times within the same case, the drug–adverse-event combination was counted only once. This analysis also included drugs recorded with drug roles such as primary suspect (PS), secondary suspect (SS), concomitant (C), or interacting (I) ( Scheme 1 ).\nUsing the PubChemPy library ( https://github.com/mcs07/pubchempy ; accessed on 13 June 2025), we searched PubChem to obtain SMILES strings for the 9400 drug names recorded in the FAERS DRUG table [ 55 , 56 ]. The constructed GnRHR agonist prediction model was then applied to 5523 drugs with SMILES strings obtained. To ensure consistency and reproducibility in descriptor calculation for numerous drugs, we limited the analysis to drugs whose SMILES structures could be directly obtained from PubChem. Drugs requiring conversion from other structural formats to SMILES were excluded. Model output probabilities were standardized relative to the cutoff value of the final prediction model using Equation 6. Structural similarity between compounds in the Tox21 database and FAERS-listed drugs predicted to have GnRHR agonist activity was evaluated visually using principal component analysis (PCA). For PCA, the 1761 molecular descriptors employed for prediction model construction were used after preprocessing. Tox21 compounds underwent PCA, and FAERS-listed drugs were projected into the same principal component space using the resulting component loadings. We considered that FAERS-listed drugs located within the range of Tox21 compounds in the PC1–PC2 plot were within the model’s chemical applicability domain (AD). As a supplementary validation of the AD, a distance-based AD assessment was additionally conducted using a k-nearest neighbors (kNN) approach. For each compound, the mean distance to its five nearest neighbors (k = 5) in the training set was calculated using the kneighbors function. The AD threshold was defined as the 95th percentile of the training-set distance distribution, corresponding to retention of 95% of the training compounds. Compounds with mean distances below the threshold were considered inside the AD, whereas those exceeding the threshold were considered outside the AD. Furthermore, the distribution of normalized scores was used as a complementary criterion for AD assessment. Several uncertainty intervals were evaluated ( Supplementary Table S4 ), and compounds with normalized active-class scores between 0.3 and 0.7 were considered outside the AD because predictions near the classification boundary were associated with greater uncertainty. The 0.3–0.7 interval was selected as a practical compromise between retaining compound coverage (86.7%) and excluding uncertain predictions. Although a normalized score of 0.5 corresponds to the model decision threshold, only compounds with normalized scores ≥0.7 were considered high-confidence positive predictions and were included in subsequent adverse event analyses. Additionally, labels were corrected according to experimental activity for 50 drugs with Tox21 GnRHR assay data whose model prediction was inconsistent with the observed experimental activity.\nWe constructed 2 × 2 contingency tables to identify adverse-event signals related to GnRHR agonists ( Table 3 ). The unit of the 2 × 2 table was the number of cases. For each MedDRA PT, cases with at least one drug judged to be positive for GnRHR agonist activity were considered to have GnRHR-active drug exposure, whereas cases involving other drugs were considered to have other drug exposure. Each case was counted only once according to the presence/absence of the adverse event of interest. Specifically,  a  was defined as the number of cases with GnRHR-active drug exposure and the adverse event of interest;  b  as the number of cases with GnRHR-active drug exposure but without the adverse event;  c  as the number of cases with other drug exposure and the adverse event; and  d  as the number of cases with other drug exposure but without the adverse event. When any cell in the 2 × 2 table contains zero, the conventional reporting odds ratio (ROR) cannot be calculated. Additionally, estimates become unstable when cell counts are very small. This bias was corrected using the Haldane–Anscombe correction, which adds 0.5 to all cells, to calculate the adjusted RORs [ 57 ]. The  p  values were calculated using Fisher’s exact test. Multiple testing was controlled using the Benjamini–Hochberg procedure. Relative reporting tendencies were assessed by detecting adverse-event signals expressed as RORs. This disproportionality analysis identified frequently reported adverse events [ 58 ]. For numerical handling, the minimum  p -value was set to 5 × 10 −324  and the minimum ROR to 1 × 10 −308 . Volcano plots were generated, with −log( q ) on the y-axis and ln(ROR) on the x-axis, for the visual interpretation of adverse-event signals.\nStatistical data, including RORs,  p  values,  q  values, and 95% confidence intervals, were analyzed using Python version 3.12.7 (Python Software Foundation, Wilmington, DE, USA). Volcano plots were generated using JMP Student Edition 19.1.2 (SAS Institute Inc., Cary, NC, USA). A random seed of 42 was used for all procedures involving stochastic processes to ensure reproducibility. The major Python packages used in this study included RDKit (version 2026.3.2), mordred-community (version 2.0.7), scikit-learn (version 1.6.1), Optuna (version 4.4.0), LightGBM (version 4.6.0), XGBoost (version 3.2.0), imbalanced-learn (version 0.13.0) and SHAP (version 0.49.1). A  q -value below 0.05 was considered statistically significant.\n\nA GnRHR agonist prediction model based on the Tox21 database was constructed in this study, achieving an ROC-AUC of 0.810 in nested cross-validation. This model was applied to FAERS-listed drugs, and associations between adverse events and drugs predicted to be GnRHR agonists were comprehensively analyzed. Results indicated that infections and infestations, respiratory, thoracic, and mediastinal disorders, and hepatobiliary disorders were more likely to be reported than other adverse events, suggesting that GnRHR agonist activity may be associated with these adverse-event reporting profiles. These findings highlight the need for early screening of GnRHR agonist activity in drug development risk assessment and postmarketing safety surveillance and may offer guidance in predicting and preventing adverse effects related to GnRHR agonist activity.","source_license":"CC-BY-4.0","license_restricted":false}