Investigation of Adverse Events Associated with Predicted GnRHR Agonists Using the FAERS Database

1. Introduction

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

GnRHR 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.

The 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].

Large-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].

In 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.

2. Results

2.1. Model Comparison

The 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. In addition, the final model retrained using all data was employed in the analysis.

Feature 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.

2.2. FAERS-Derived Drug Screening

In total, 5,523 FAERS-derived 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, 5,177 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, 1,191 drugs remained within the applicability domain, of which 367 were predicted to be positive for GnRHR agonist activity (Supplementary Table 1).

To 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, 8,066 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 agonist–positive drugs 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.

2.3. Adverse-Event Signals

First, 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 531 important adverse events met the following criteria: at least 1,000 reports, ln(ROR) >1, and p < 0.05 (Supplementary Table 2). When adverse events with at least 1,000 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 1), infections and infestations (Supplementary Figure 2), and hepatobiliary disorders (Supplementary Figure 3).

3. Discussion

3.1. Applicability of the GnRHR Agonist Prediction Model

Using GnRHR agonist activity data from the Tox21 database, we developed a machine learning model based on 1,761 molecular descriptors. With an ROC-AUC of 0.810 in nested cross-validation, this model achieved reliable predictive performance. 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, and 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 [38]. 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 [39]. 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.

In 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. Therefore, GnRHR agonists may be recognized by GnRHR through 3D structures containing branched or cyclic frameworks and local electronic interactions mediated by hydrogen bonding rather than primarily through hydrophobic or pi interactions. Indeed, GnRHR contains many polar amino acids and depends on hydrogen bonding and electrostatic interactions [40]. 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 [41,42]. 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 align 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.

Drugs 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-mediated adverse event risks may be evaluable not only for established GnRHR agonists but also for drugs with potential GnRHR agonist activity.

3.2. Adverse Events Related to Respiratory, Thoracic, and Mediastinal Disorders

Drugs 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.

Several cohort and database studies have reported that GnRH agonist therapy is associated with increased risks of interstitial lung disease [11], community-acquired pneumonia [43], and radiation pneumonitis [44].

Ojeda et al. assessed morphological and biochemical changes in the lungs of male rats receiving surgical or pharmacological castration and found that androgen deficiency induced alveolar septal destruction, emphysema-like changes, and alterations in the phospholipid composition of pulmonary surfactant [45]. 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.

Taken together, androgen deficiency induced by GnRH agonists may compromise the structural integrity of lung tissue through morphological and biochemical changes in pulmonary cell membranes, thereby increasing susceptibility to respiratory diseases, including pneumonia.

3.3. Adverse Events Related to Infections and Infestations

The high number of reported infection-related adverse events aligns with a series of evidence linking androgens to immune function.

Androgen receptors are expressed in many immune cells, including neutrophils [46], and androgens promote granulocyte progenitor proliferation and increase neutrophil production [47]. Using androgen receptor–knockout mice, Chuang et al. showed that androgen receptor deficiency was associated with neutropenia and increased vulnerability to microbial infection and that exogenous androgen administration restored neutrophil counts in castrated mice [48].

According to Ojeda et al., androgen deficiency induces changes in the phospholipid composition of pulmonary surfactant and morphological changes in lung tissue [45], 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 [43]. Thus, the synergistic effects of immune dysfunction and pulmonary tissue changes may increase infection risk among GnRH agonist users.

3.4. Adverse Events Related to Hepatobiliary Disorders

Hepatobiliary adverse events were also frequently reported for drugs predicted to be GnRH agonists in this study, consistent with the known metabolic adverse effects of ADT.

In a large population-based study, Gild et al. showed that ADT, including GnRH 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 GnRH agonist therapy significantly increased biliary disease incidence and that the risk increased with longer duration of use, suggesting a dose–response relationship [49]. In addition, bilateral orchiectomy was not associated with a significant increase in biliary disease risk, suggesting that mechanisms specific to GnRH agonists contribute to biliary disease.

Moreover, GnRH agonists induce metabolic abnormalities, including increased fat mass, particularly abdominal fat accumulation, elevated serum triglycerides, and increased insulin resistance [50]. According to a metabolomic analysis, most bile acids and their metabolites increase within 3 months after the initiation of GnRH agonist therapy [51]. Additionally, hormonal changes can induce gallbladder dysmotility [52]. These factors may have contributed to the increased risk of developing hepatobiliary diseases observed in the present study.

3.5. Adverse Events Related to the Reproductive System and Breast Disorders

In the present study, the proportion of reported adverse events classified as reproductive system and breast disorders was relatively low in the GnRHR agonist group (Supplementary Figure 4). 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], changes in the number of reports related to these diseases may have influenced the overall reporting structure, making the proportion of reported reproductive system and breast disorders appear relatively low.

Conversely, ovarian hyperstimulation syndrome (ROR = 13.48, p < 5 × 10⁻³²⁴) and polycystic ovaries (ROR = 3.17, p = 1.76 × 10^-79), 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.

Therefore, rather than being interpreted as an indication of a reduced risk of reproductive system–related adverse events, the lower proportion of reported reproductive system and breast disorders may reflect the influence of underlying indications and changes in the reporting structure.

3.6. Limitations

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

Second, 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. This limitation affects how the prediction results are interpreted clinically. Third, the model was applied to 1,191 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.

Fourth, 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 [54,55]. 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.

Fifth, spontaneous reporting databases record only cases with reported adverse events and without information on all patients exposed to a drug [56]. 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 [56,57].

Finally, when multiple drugs are used concomitantly, identifying which drug caused an adverse event is difficult [56,58]. 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.

4. Materials and Methods

4.1. Data Source

We collected data on GnRHR agonists from the Tox21 database (AID: 1920067; https://pubchem.ncbi.nlm.nih.gov/bioassay/1920067; accessed January 6, 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, 7,133 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, including inconclusive compounds, were considered negative (label: 0). Ultimately, the final dataset comprised 7,071 negative and 62 positive compounds.

4.2. Molecular Structures and Chemical Descriptors

RDKit, an open-source cheminformatics toolkit, was used for cleaning and standardizing SMILES strings [18] 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 [19], followed by energy minimization with MMFF [20]. The selected conformer was then converted into structure-data file (SDF) format.

For each compound, molecular descriptors were calculated using the Python libraries RDKit and Mordred [21,22]. Two-dimensional and 3D descriptors, including RDKit descriptors (221 dimensions) and Mordred descriptors (1,826 dimensions), were obtained. From a total of 2,047 descriptors, perfectly correlated descriptors and those with zero variance were removed, leaving 1,761 descriptors for model construction.

4.3. Model Development

Classification models were developed to predict GnRHR agonist activity. Four algorithms were compared: RandomForest [23], LightGBM [24], XGBoost [25], and BalancedRandomForest [26]. The 1,761 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.

As 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 [27]. 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 [28].

4.4. Evaluation Metrics

The 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.

Sensitivity (SE): SE is the proportion of actual positive samples correctly predicted as positive.

$$SE={\displaystyle \frac{TP}{TP+FN}}$$

(1)

Specificity (SP): SP is the proportion of actual negative samples correctly predicted as negative.

$$SP={\displaystyle \frac{TN}{TN+FP}}$$

(2)

Accuracy (ACC): ACC is the number of correctly predicted samples divided by the total number of samples.

$$ACC={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(3)

Balanced accuracy (BA): BA is the mean of SE and SP.

$$BA={\displaystyle \frac{1}{2}}\left(SE+SP\right)$$

(4)

Matthews correlation coefficient (MCC): MCC is a metric for evaluating classification performance in imbalanced datasets [29].

$$MCC={\displaystyle \frac{TP\times TN-FP\times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}}$$

(5)

AUC: AUC is calculated as the area under the ROC curve based on SE and 1 - specificity (1 - SP) across thresholds [30].

Average precision (AP): AP is calculated as the area under the precision-recall curve based on precision [TP/(TP + FP)] and recall (SE) across thresholds [31].

The optimal cutoff value for TP, FN, TN, and FP was determined using Youden’s index, which maximizes the sum of SE and SP [32]. Given that the final prediction model had a model-specific cutoff, the model output probabilities were standardized for drug prediction using the following equation:

$$X_n=X_u^{-{\log }_{c}2}$$

(6)

where 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. This transformation standardizes the model output probabilities relative to the cutoff value, allowing for consistent interpretation of prediction scores for FAERS-listed drugs.

4.5. FAERS Database

The FAERS database was accessed from the official FDA website [33], 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. Overview of the FAERS data integration and disproportionality analysis workflow. SMILES, Simplified Molecular Input Line Entry System; ID, identification; GnRHR, gonadotropin-releasing hormone receptor; ln(ROR), natural logarithm of the reporting odds ratio.

Scheme 1. Overview of the FAERS data integration and disproportionality analysis workflow. SMILES, Simplified Molecular Input Line Entry System; ID, identification; GnRHR, gonadotropin-releasing hormone receptor; ln(ROR), natural logarithm of the reporting odds ratio.

4.6. Application of the Prediction Model to FAERS-listed Drugs

Using the PubChemPy library (https://github.com/mcs07/pubchempy; accessed June 13, 2025), we searched PubChem to obtain SMILES strings for the 9,400 drug names recorded in the FAERS DRUG table [34,35]. The constructed GnRHR agonist prediction model was then applied to 5,523 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; thus, 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. Drugs with normalized scores of 0.7 or higher were classified as predicted positive for GnRHR agonist activity. 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 1,761 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. 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.

4.7. Association Analysis Between Predicted GnRHR Agonist Activity and Adverse Events

We 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 [36]. The p values were calculated using Fisher’s exact test. Relative reporting tendencies were assessed by detecting adverse-event signals expressed as RORs. This disproportionality analysis identified frequently reported adverse events [37]. For numerical handling, the minimum p-value was set to 5 × 10⁻³²⁴ and the minimum ROR to 1 × 10⁻³⁰⁸. Volcano plots were generated, with -log(p) on the y-axis and ln(ROR) on the x-axis, for the visual interpretation of adverse-event signals.

4.8. Statistical Analysis

Statistical data, including RORs, p values, and 95% confidence intervals, were analyzed using Python scripts (Python version 3.12.7). Volcano plots were generated using JMP Student Edition 19.1.2 (SAS Institute Inc., Cary, NC, USA). A p-value below 0.05 was considered statistically significant.

5. Conclusions

A GnRHR agonist prediction model based on the Tox21 database was constructed in this study, achieving an 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.