# Methods — Data-Driven Reclassification of Pharmacogenomic ADR Associations *Living document. Describes the current accepted model only (no development history). Update after each accepted change to features, clustering, or labels.* ## Overview We reclassified gene × drug × ICD10 adverse-drug-reaction (ADR) associations (207 locus × phenotype signals after collapsing co-located variants; see below) from a rule-based four-tier expert cascade into data-driven clusters. The objective was to test whether the expert tiers (Indication, Comorbidity, Predisposition, Effect) correspond to natural structure in the molecular evidence, and to expose borderline associations that a hard-threshold cascade forces into a single tier. The reclassification uses **only raw measurements and continuous or named statistical transformations of them** — no threshold or output of the expert cascade is used as an input, so agreement with the cascade is an unbiased diagnostic rather than a built-in result. ## Association unit (locus collapse) The unit of analysis is one **locus × phenotype** signal, not one variant. Within each (drug, ICD10) phenotype, co-located variants (same chromosome, within a 250 kb window) are collapsed to a single representative — the lead variant with the strongest drug-ADR association (`−log10 P`) — because variants in linkage disequilibrium at one locus are not independent associations and would otherwise inflate cluster sizes and apparent replication. Gene-level burden tests (no single position) are each retained as their own locus. This collapses the raw variant-level table to 207 locus × phenotype signals; the number of variants merged per locus and their genes are retained as reported annotations. ## Feature construction Each association is represented by 27 features. All features are continuous (or ordinal integer) and standardized to zero mean and unit variance before clustering. **No feature is weighted relative to any other**: every feature enters the distance metric at unit variance, so the clustering reflects the data rather than analyst preference. The guiding principle of feature selection was *parsimony with justification*: a feature is retained only if it contributes non-redundant structure, established by leave-one-out cluster-stability testing and by recovery of the independent literature-anchor benchmark — not by overall agreement with the expert cascade. Features that were noisy copies of a retained feature, or that re-introduced estimation instability, were removed; features that reinforce a real evidence axis as an ensemble were kept even when mutually correlated. ### Genome-wide association backbone For every variant we extracted association statistics from four GWAS phenotypes run on the UK Biobank cohort, providing four orthogonal views of the same locus: | phenotype | tool | question it answers | |-----------|------|---------------------| | drug-ADR | REGENIE (discovery) | does the variant associate with the adverse event under drug exposure? | | disease (ICD10) | REGENIE | does it associate with the disease in the general population, independent of the drug? | | drug dose | REGENIE | does it affect the dose patients are titrated to? | | drug prescribed | REGENIE | does it affect who is prescribed the drug (a proxy for the indication)? | For each phenotype we use the absolute effect size (`log1p|BETA|`) and the significance (`−log10 P`). The standard error is retained only for the discovery (drug-ADR) arm; the standard errors of the other three arms were a single redundant axis dominated by rare-variant estimation instability (Firth-Wald quasi-separation in burden tests) and were removed. Standard errors and odds ratios are winsorized at the 95th percentile for the same reason — to prevent a small number of quasi-separated burden tests from dominating the Euclidean distance metric. ### Variant properties A single feature, the log carrier count in cases (`log1p(A1_CASE_CT)`), captures both allele rarity (carrier count scales with minor-allele frequency) and statistical power: associations with zero carriers in cases (≈23% of variants) have no power to be confirmed, a signal that allele frequency alone cannot convey. Minor-allele frequency was tested as a separate feature and removed as a noisy copy of the carrier count. ### Independent-replication evidence A held-out PLINK2 validation analysis provides an independent replication of the discovery signal. We summarize it by the validation `−log10 P`. In a discovery→validation design the discovery arm is selected for significance and is therefore always the stronger arm; the validation p-value is consequently the *weakest link* (minimum statistic) of the two arms, which is the appropriate conjunction measure of confirmation: an association is confirmed only insofar as its weaker arm remains significant. A difference between the two arms' p-values was rejected as a replication measure because it inverts the desired behaviour (a large but still-significant drop registers as weak replication). Effect- direction consistency between the discovery BETA and validation odds ratio is reported as a supplementary column but excluded from clustering, as it describes an axis (rare-variant winner's-curse shrinkage) orthogonal to the classification. ### Cross-phenotype effect contrasts Four signed contrasts compare the drug-ADR effect against the other phenotypes, encoding drug-specific amplification and indication confounding without any threshold (in contrast to the cascade's ratio cutoffs): discovery-vs-disease, discovery-vs-dose, discovery-vs-prescribed, and disease-vs-prescribed. The three contrasts among the binary (logistic) arms — discovery, disease and prescribed — are differences of `log1p|BETA|` (log-odds ratios), which are comparable across arms up to a constant scale offset that global z-scoring removes. The drug-dose phenotype, however, is a **linear** trait (ln median dose), so its coefficient is not on the log-odds scale; the discovery-vs-dose contrast is therefore computed from **within-arm-standardized** per-allele effects (each arm's allele-frequency- scaled effect `BETA·√(2f(1−f))` z-scored across variants), making the two arms unit-free and comparable. This avoids the scale artifact whereby the rare-variant discovery arm's inflated effect sizes (median |BETA| ≈ 4.6 vs ≈ 0.02 for the linear dose arm) would otherwise dominate the contrast. Because the contrasts above use absolute effect sizes, a separate **effect-direction** feature retains the sign: the product of the discovery and disease effect signs (+1 when the variant pushes the disease the same way with and without the drug — vertical pleiotropy consistent with a real mechanism; −1 when the drug-arm effect opposes the general-population effect, the signature of a protective or confounded association). This distinguishes drug-caused from drug-protective associations, which absolute-magnitude features cannot. ### Pharmacovigilance and internal pleiotropy Real-world adverse-event signal is captured by FDA Adverse Event Reporting System (FAERS) disproportionality metrics — the Bayesian Information Component lower bound and the log Proportional Reporting Ratio. Variant pleiotropy is captured by an **internal pleiotropy** count: the number of the project's own GWAS phenotypes in which the variant reaches significance. This is computed entirely from the GWAS corpus (no external database) and is variant-level rather than gene-level: a variant significant only in the drug-ADR phenotype is specific (a candidate), whereas one significant across many phenotypes is non-specific. Effect-tier variants carry the lowest internal pleiotropy, consistent with drug-specific signals. ### Self-contained design (no external gene annotations) The model deliberately uses no external, hand-curated, or queried gene annotations: tissue expression (GTEx/HPA), druggability (DGIdb), external PheWAS/GWAS-Catalog pleiotropy, and gene-disease association scores were evaluated but excluded from the feature set, so the classification is fully reproducible from the GWAS corpus and generalises genome-wide. Gene-disease plausibility (Open Targets) is retained only as a reported annotation for the candidate shortlist, not as a clustering input. ### Combined-evidence statistics Two orthogonal evidence layers — pharmacovigilance (FAERS) and statistical replication — are summarized by two named meta-analytic statistics, each z-scored: - **Stouffer's combined Z** (Σz/√k; Stouffer, 1949) — the overall evidence strength across layers, high when any layer is strong. - **Minimum-statistic for conjunction inference** (min of the two z-scores; Nichols et al., 2005) — high only when both layers are jointly elevated. This is the continuous, threshold-free analogue of the cascade's count of independent evidence layers, and is the single feature that positively marks the Effect class. One further interaction completes the set: a continuous drug-specificity log-ratio (discovery vs. disease effect), the threshold-free analogue of the amplification ratio. ### Temporal indication signal The temporal relationship between prescription and diagnosis is captured by the fraction of ICD10 events occurring after first prescription and by a directional log-ratio of post- to pre-prescription event fractions. A low ratio (disease preceded the drug) sharply marks the Indication class. Membership of the ICD10 code on the drug's official indication list is included as a categorical feature and is the single strongest Indication marker. A binary flag marks non-specific outcomes (ICD10 R-chapter symptom codes and Z-chapter health-status codes), which are of low interpretability and account for a disproportionate share of the nominal Effect tier. ## Clustering The standardized feature matrix was clustered with four methods spanning model families — k-means, Gaussian mixture (tied covariance), latent class analysis, and HDBSCAN — with the number of clusters selected per method by silhouette (k-means), Bayesian Information Criterion (mixture, latent class), or density (HDBSCAN). k-means and latent class analysis are reported as the primary methods (both fully stable under consensus), with the mixture and density methods as robustness checks; on the unweighted feature set the four methods agree closely (≈55% concordance with the expert tiers at four clusters, with consistent per-class recovery). For the stochastic methods, 100-seed consensus clustering was used to derive soft membership probabilities and a per-association stability score (the maximum consensus membership), so that borderline associations are reported as borderline rather than forced into a single class — directly addressing the reviewer requirement for a differentiating rather than segregating classifier. ## Class labelling (EPIC) Clusters are labelled against the four-class EPIC scheme using the 18-pair benchmark of literature-validated positive (ADR) and negative (indication) associations as anchors, by majority vote of anchor membership: - **E — Effect**: full mechanism; all evidence layers jointly elevated. - **P — Predisposition**: drug-amplified, often rare-variant signal. - **I — Indication**: on-label association; disease precedes the drug. - **C — inConclusive**: the low-evidence residual. The data show this class is not a positive comorbidity signature — it is negative on every evidence axis, with only ~28% of members showing a genuine comorbidity pattern (independent disease association without drug amplification) and the majority weak across all layers. We therefore relabel it *inconclusive*, the category of associations that do not clear the evidence threshold for any positive class. ## Scope of the unsupervised classifier and the Effect class The unsupervised methods are a robust discriminator of **Indication**, **inConclusive** and **Predisposition** structure, but they do **not** recover the **Effect** class: on a curated gold standard of 35 pharmacologically-defined pairs, they achieved 0.15–0.26 accuracy and zero Effect recall, collapsing clinically-distinct associations onto the dominant rare-variant Predisposition axis. The expert rule cascade, by contrast, recovers Indication with F1 0.96 and inConclusive with F1 0.82, but is itself conservative on Effect (recall 0.22 — genuine adverse reactions such as ACE-inhibitor angioedema and hyponatraemia, gabapentin dizziness, and antiplatelet haemorrhage are frequently assigned to Predisposition). We therefore do not present the unsupervised clustering as a stand-alone replacement for the four-class scheme. **Effect identification requires expert or biological input**: it is flagged jointly by the rule cascade, by the effect-direction (sign-concordance) feature — which distinguishes drug-amplified vertical pleiotropy from protective or confounded associations — and by an Open-Targets gene-disease plausibility score reported per association (the Effect tier carries the lowest mean gene-disease plausibility, 0.002, confirming its gene-level signals are largely incidental). The data-driven clustering is used to expose borderline associations and to reproduce the Indication / inConclusive / Predisposition structure without the rule cascade's hard thresholds, not to adjudicate Effect. ## Reported (non-clustering) annotations Two annotations are computed and reported per association but excluded from the feature matrix, to avoid circularity (external knowledge) or because they describe an axis orthogonal to the classification: the Open-Targets gene-disease plausibility score, and a per-locus variant-multiplicity count flagging LD-clustered signals that would otherwise inflate apparent independent replication. ## Software and data Outputs are written additively (new columns and sheets) to the atlas workbook; the source classification is never modified. GWAS summary statistics are from the UK Biobank PGxPred analysis; tissue expression from GTEx v8 and the Human Protein Atlas; ICD10→tissue mapping from manual curation with MONDO/UBERON fallback.