People who carry rare loss-of-function variants in a gene are a natural experiment: they show what happens when the protein that gene encodes stops working for a lifetime. When that protein is also a drug target, their phenotypes are a preview of what the drug might do. Drug discovery teams use this idea to support targets and to anticipate on-target adverse effects before clinical data exist.
The Human Phenotype Ontology (HPO)
annotates genes with the clinical features of the Mendelian diseases
they cause. dbparser reads the HPO gene-phenotype
associations distributed by Harmonizome
with parseHPO(). It then links them to DrugBank targets by
gene symbol with merge_drugbank_hpo().
This vignette asks one question for three drugs:
Do the Mendelian phenotypes of a drug’s target genes agree with what we already know about the drug’s therapeutic and adverse effects?
The three drugs act through different mechanisms:
| Drug | Main target | Use |
|---|---|---|
| Bivalirudin | Thrombin (F2) | Anticoagulant during coronary interventions |
| Leuprolide | GnRH receptor (GNRHR) | Androgen deprivation, endometriosis |
| Calcitriol | Vitamin D receptor (VDR) | Hypocalcaemia, renal bone disease |
The package ships a subset of the Harmonizome HPO file with every
association for seven genes. To analyse all genes, download
gene_attribute_edges.txt.gz from Harmonizome and pass its
path to parseHPO(); gzip-compressed files are read
directly.
hpo_db <- parseHPO(
system.file("extdata", "hpo_gene_attribute_edges.txt", package = "dbparser"),
db_version = "Harmonizome 2015")
hpo_db$genes
#> # A tibble: 7 × 3
#> gene_symbol entrez_id n_phenotypes
#> <chr> <chr> <int>
#> 1 F2 2147 71
#> 2 GNRH1 2796 24
#> 3 GNRHR 2798 30
#> 4 KCNQ2 3785 69
#> 5 MPO 4353 7
#> 6 SCN1A 6323 101
#> 7 VDR 7421 183search_phenotypes() and
get_genes_for_phenotype() help explore the phenotypes:
search_phenotypes(hpo_db, "bleeding")
#> # A tibble: 3 × 4
#> hpo_id phenotype_name n_genes genes
#> <chr> <chr> <int> <chr>
#> 1 HP:0000225 gingival bleeding 1 F2
#> 2 HP:0001892 abnormal bleeding 1 F2
#> 3 HP:0003010 prolonged bleeding time 1 F2
get_genes_for_phenotype(hpo_db, "hypogonadotrophic hypogonadism", exact = TRUE)
#> # A tibble: 2 × 5
#> gene_symbol entrez_id hpo_id phenotype_name association_weight
#> <chr> <chr> <chr> <chr> <dbl>
#> 1 GNRH1 2796 HP:0000044 hypogonadotrophic hypogon… 1
#> 2 GNRHR 2798 HP:0000044 hypogonadotrophic hypogon… 1Leuprolide and calcitriol come from a saved DrugBank sample. Bivalirudin is parsed from a DrugBank XML record.
two_drugs <- readRDS(system.file("two_drugs.RDS", package = "dbparser"))
bivalirudin <- parseDrugBank(
db_path = system.file("extdata",
"drugbank_record_small_molecule.xml",
package = "dbparser"),
drug_options = NULL,
parse_salts = FALSE,
parse_products = FALSE,
references_options = NULL,
cett_options = "targets")merge_drugbank_hpo() maps each DrugBank target to its
gene and adds drugbank_id and drug_name to the
HPO associations. We merge each DrugBank object and stack the
results.
merged <- lapply(list(two_drugs, bivalirudin), merge_drugbank_hpo,
hpo_db = hpo_db)
# drug -> target -> gene bridge, including targets without HPO data
gene_bridge <- bind_rows(lapply(merged, function(m) {
m$integrated_data$DrugBank_GeneSymbol_Mapping
}))
# HPO associations of the drug target genes
target_phenotypes <- bind_rows(lapply(merged, function(m) {
m$integrated_data$gene_phenotype_enriched
})) %>%
filter(!is.na(drugbank_id)) %>%
distinct(drug_name, gene_symbol, hpo_id, phenotype_name)
gene_bridge %>%
select(drug_name, target_id, polypeptide_name, gene_name) %>%
left_join(count(target_phenotypes, gene_symbol, name = "hpo_terms"),
by = c("gene_name" = "gene_symbol"))
#> # A tibble: 5 × 5
#> drug_name target_id polypeptide_name gene_name hpo_terms
#> <chr> <chr> <chr> <chr> <int>
#> 1 Leuprolide BE0000203 Gonadotropin-releasing hormone rece… GNRHR 30
#> 2 Leuprolide BE0010526 Progonadoliberin-1 GNRH1 24
#> 3 Calcitriol BE0000779 Vitamin D3 receptor VDR 183
#> 4 Calcitriol BE0008987 Homeobox protein Hox-A10 HOXA10 NA
#> 5 Bivalirudin BE0000048 Prothrombin F2 71Four of the five targets have HPO annotations. HOXA10, a
target that DrugBank lists for calcitriol with an unknown role
(known_action = "unknown"), has none. Missing phenotypes
are a coverage gap, not evidence that the target is safe to
modulate.
GNRH1 encodes the GnRH precursor, the natural ligand of the GnRH receptor. Leuprolide is a GnRH analogue, so GNRH1 and GNRHR describe the same signalling axis.
Harmonizome propagates every annotation up the ontology. A gene annotated with epistaxis is also annotated with abnormality of the nose, phenotypic abnormality and the root term all. These grouping terms, the inheritance and onset terms, and the clinical modifiers say nothing about the gene, so we drop them.
is_uninformative <- function(hpo_id, phenotype_name) {
hpo_id %in% c("HP:0000001", "HP:0000118") |
grepl("inheritance|onset|clinical modifier|variability|expressivity",
phenotype_name) |
grepl("^(functional )?abnormal(ity)? ", phenotype_name) |
grepl(" abnormality$|anomaly", phenotype_name)
}
specific_phenotypes <- target_phenotypes %>%
filter(!is_uninformative(hpo_id, phenotype_name))
count(target_phenotypes, drug_name, gene_symbol, name = "all_terms") %>%
left_join(count(specific_phenotypes, drug_name, gene_symbol,
name = "specific_terms"),
by = c("drug_name", "gene_symbol"))
#> # A tibble: 4 × 4
#> drug_name gene_symbol all_terms specific_terms
#> <chr> <chr> <int> <int>
#> 1 Bivalirudin F2 71 21
#> 2 Calcitriol VDR 183 70
#> 3 Leuprolide GNRH1 24 6
#> 4 Leuprolide GNRHR 30 11This name-based filter is a quick heuristic. For a full analysis,
filter on the ontology structure instead, e.g. by keeping only the most
specific terms with the ontologyIndex package and the HPO
.obo file.
To compare targets, we assign each phenotype to a clinical domain with keyword rules. The rules are applied in order, and the first match wins.
domain_rules <- c(
"Bleeding & coagulation" = "hemorrhag|bleed|epistaxis|menorrhagia|bruising|prothrombin|thromboplastin",
"Thrombosis" = "thrombo(sis|phlebitis)|embolism",
"Reproductive & gonadal" = "gonad|puberty|genital|test[ei]|amenorrh|menstrua|androgen",
"Calcium & mineral" = "calc|phosphat|parathyroid|vitamin|alkaline|nephrolith",
"Skeletal & dental" = "bone|skelet|rickets|bowing|bossing|fractur|osteo|epiphys|metaphys|fontanel|sutur|rib|costochondral|genu|scoliosis|teeth|dental|enamel|occiput|cephaly|wrist|ankle|trabecul|cortex|stature",
"Neuromuscular & growth" = "seizure|hypotonia|motor|gait|muscle|walking|standing|irritab|neurodevelop|anosmia|growth|thrive|body weight"
)
assign_domain <- function(phenotype_name) {
domain <- rep("Other", length(phenotype_name))
for (rule in rev(names(domain_rules))) {
domain[grepl(domain_rules[[rule]], phenotype_name)] <- rule
}
domain
}
specific_phenotypes <- specific_phenotypes %>%
mutate(domain = assign_domain(phenotype_name),
drug_target = paste0(drug_name, " (", gene_symbol, ")"))
domain_counts <- specific_phenotypes %>%
count(drug_target, domain)
domain_counts
#> # A tibble: 11 × 3
#> drug_target domain n
#> <chr> <chr> <int>
#> 1 Bivalirudin (F2) Bleeding & coagulation 12
#> 2 Bivalirudin (F2) Other 1
#> 3 Bivalirudin (F2) Reproductive & gonadal 2
#> 4 Bivalirudin (F2) Thrombosis 6
#> 5 Calcitriol (VDR) Calcium & mineral 10
#> 6 Calcitriol (VDR) Neuromuscular & growth 15
#> 7 Calcitriol (VDR) Other 5
#> 8 Calcitriol (VDR) Skeletal & dental 40
#> 9 Leuprolide (GNRH1) Neuromuscular & growth 2
#> 10 Leuprolide (GNRH1) Reproductive & gonadal 4
#> 11 Leuprolide (GNRHR) Reproductive & gonadal 11The plot below shows each target’s phenotype profile.
domain_levels <- c(names(domain_rules), "Other")
tidy_to_matrix <- function(counts) {
targets <- sort(unique(counts$drug_target))
m <- matrix(0, nrow = length(domain_levels), ncol = length(targets),
dimnames = list(domain_levels, targets))
m[cbind(counts$domain, counts$drug_target)] <- counts$n
as.data.frame(m)
}
canvasXpress(
data = tidy_to_matrix(domain_counts),
graphType = "Stacked",
graphOrientation = "horizontal",
title = "Mendelian phenotypes of drug target genes",
subtitle = "Specific HPO terms per target, by clinical domain",
xAxisTitle = "Number of HPO terms",
colorScheme = "Tableau",
legendPosition = "right",
showSampleNames = TRUE
)Each target has a distinct profile: F2 is dominated by bleeding and thrombosis, the GnRH axis genes by reproductive phenotypes, and VDR by skeletal and mineral phenotypes.
HPO tells us what happens when a gene is altered, mostly by loss of function. It does not tell us whether a drug mimics or opposes that change. That depends on the drug’s action on the target.
For bivalirudin, the action comes from the parsed DrugBank record.
Each CETT action row carries its drugbank_id, so we can
read the action for this drug alone:
bivalirudin$cett$targets$actions %>%
left_join(select(bivalirudin$cett$targets$general_information,
target_id, drugbank_id, target_name = name),
by = c("target_id", "drugbank_id"))
#> # A tibble: 1 × 4
#> action target_id drugbank_id target_name
#> <chr> <chr> <chr> <chr>
#> 1 inhibitor BE0000048 DB00006 ProthrombinThe two_drugs sample was saved before CETT actions
carried a drugbank_id, so its actions table cannot be
attributed to a single drug. For leuprolide and calcitriol we write down
the mode of action from their DrugBank entries and their known
pharmacology:
drug_direction <- tibble::tribble(
~drug_name, ~gene_symbol, ~action, ~net_effect_on_target,
"Bivalirudin", "F2", "inhibitor", "loss of function",
"Leuprolide", "GNRHR", "agonist", "loss of function (desensitisation)",
"Leuprolide", "GNRH1", "analogue", "loss of function (desensitisation)",
"Calcitriol", "VDR", "agonist", "gain of function"
)
drug_direction
#> # A tibble: 4 × 4
#> drug_name gene_symbol action net_effect_on_target
#> <chr> <chr> <chr> <chr>
#> 1 Bivalirudin F2 inhibitor loss of function
#> 2 Leuprolide GNRHR agonist loss of function (desensitisation)
#> 3 Leuprolide GNRH1 analogue loss of function (desensitisation)
#> 4 Calcitriol VDR agonist gain of functionLeuprolide is a receptor agonist, but continuous exposure desensitises the GnRH receptor and shuts down the pituitary–gonadal axis. Its net effect is therefore similar to losing the receptor.
Now we look at the domains that best represent each target:
key_domains <- c("Bleeding & coagulation", "Thrombosis",
"Reproductive & gonadal", "Calcium & mineral")
specific_phenotypes %>%
filter(domain %in% key_domains) %>%
group_by(drug_target, domain) %>%
summarise(n_terms = n(),
examples = paste(head(sort(phenotype_name), 5), collapse = "; "),
.groups = "drop") %>%
arrange(drug_target, desc(n_terms)) %>%
knitr::kable()| drug_target | domain | n_terms | examples |
|---|---|---|---|
| Bivalirudin (F2) | Bleeding & coagulation | 12 | bruising susceptibility; epistaxis; gastrointestinal hemorrhage; gingival bleeding; internal hemorrhage |
| Bivalirudin (F2) | Thrombosis | 6 | cerebral venous thrombosis; deep venous thrombosis; pulmonary embolism; recurrent thrombophlebitis; thrombophlebitis |
| Bivalirudin (F2) | Reproductive & gonadal | 2 | menstrual irregularities; puberty and gonadal disorders |
| Calcitriol (VDR) | Calcium & mineral | 10 | elevated alkaline phosphatase; elevated circulating parathyroid hormone (pth) level; hyperparathyroidism; hypocalcemia; hypocalcemic seizures |
| Leuprolide (GNRH1) | Reproductive & gonadal | 4 | delayed puberty; hypogonadism; hypogonadotrophic hypogonadism; puberty and gonadal disorders |
| Leuprolide (GNRHR) | Reproductive & gonadal | 11 | amenorrhea; androgen insufficiency; aplasia/hypoplasia of the testes; decreased testicular size; external genital hypoplasia |
The two reproductive terms for F2 (menstrual irregularities and puberty and gonadal disorders) are most likely ancestors of menorrhagia that the propagation added. They reflect bleeding, not a separate reproductive effect.
Bivalirudin inhibits thrombin, so it mimics loss of F2 function. Prothrombin deficiency causes bleeding: epistaxis, gingival and gastrointestinal haemorrhage, and prolonged prothrombin and partial thromboplastin times. Bleeding is the main adverse effect of bivalirudin, and a prolonged clotting time is how its effect is monitored. The thrombosis terms come from the gain-of-function variant (prothrombin G20210A). They describe the state that the drug is given to oppose, not a predicted adverse effect.
Loss of GNRHR or GNRH1 causes hypogonadotropic hypogonadism, androgen insufficiency, amenorrhoea and delayed puberty. Leuprolide causes the same state in patients on purpose. It is the basis of androgen deprivation therapy in prostate cancer and of its use in endometriosis and central precocious puberty. It also explains the drug’s hypogonadal adverse effects. GNRH1 alone also adds anosmia. Anosmia belongs to the Kallmann form of hypogonadotropic hypogonadism, where GnRH neurons and olfactory neurons fail to develop. The term is probably attached to GNRH1 through a disease entry that covers hypogonadism “with or without anosmia”. It reflects development, not GnRH signalling, so it is not expected from leuprolide. Checks like this are part of reading gene-level annotations.
Loss of VDR causes hereditary vitamin D–resistant rickets: rickets, hypocalcaemia, hypophosphataemia, secondary hyperparathyroidism and alopecia. Calcitriol activates the receptor, so these phenotypes do not predict its adverse effects. They describe the deficiency state that calcitriol treats. The main adverse effects of calcitriol, hypercalcaemia and hypercalciuria, come from too much receptor activity, and the loss-of-function annotations do not contain them.
The bar chart shows how many terms each target has; VDR dominates because vitamin D-resistant rickets is a multi-system disease that is described in detail. To compare the composition of the profiles, we convert the counts to the share of each target’s specific terms that fall in each domain.
domain_share <- domain_counts %>%
group_by(drug_target) %>%
mutate(n = round(100 * n / sum(n))) %>%
ungroup()
canvasXpress(
data = as.data.frame(t(tidy_to_matrix(domain_share))),
graphType = "Heatmap",
title = "Phenotype profile of each drug target",
subtitle = "Share of specific HPO terms (%)",
colorSpectrum = c("white", "#08519c"),
showHeatmapIndicator = TRUE,
heatmapIndicatorPosition = "right",
heatmapIndicatorHistogram = FALSE,
samplesClustered = FALSE,
variablesClustered = FALSE,
showSampleNames = TRUE,
showVariableNames = TRUE
)The profiles hardly overlap. Each target points to one organ system, and that system matches the drug’s known pharmacology: haemostasis for bivalirudin, the reproductive axis for leuprolide, and bone and mineral metabolism for calcitriol. The two GnRH axis genes have similar profiles, as expected for a ligand and its receptor.
Linking DrugBank targets to HPO phenotypes takes three calls:
parseHPO(), parseDrugBank() and
merge_drugbank_hpo(). In this example the Mendelian
phenotypes of the target genes agreed with known pharmacology once the
direction of drug action was taken into account:
Before using these annotations to anticipate effects for a new compound, check the direction of action. Also treat genes without annotations (HOXA10 here) as unknown, not as safe.