Explaining Drug–Drug Interactions with LINCS L1000 Signatures and Multi-Source Phenotype Gene Sets: A Workflow and Pre-Specified Case Study (Tacrolimus–Azole Nephrotoxicity)

Abstract

Background. Drug–drug interactions (DDIs) are a leading cause of adverse drug events (ADEs). For many interacting pairs, a pharmacokinetic mechanism such as CYP3A inhibition explains why exposure rises. It does not explain why a particular adverse event appears, or why only some patients are affected. Text-mining studies have catalogued thousands of DDI pairs together with the ADEs reported for them¹, but most of these pairs have no gene-level explanation.

Approach. We describe a reusable workflow, "Explain Drug–Drug Interactions (LINCS)", that proposes gene-level explanations for a DDI and its adverse event. It integrates two kinds of gene sets. The first comes from drug-induced transcriptional perturbations: the LINCS L1000 consensus up- and down-regulated genes for each drug². The second comes from the adverse-event phenotype: human and mouse phenotype annotations (HPO, MGI-MP), literature gene mentions (Geneshot³) and GWAS Catalog associations. The workflow intersects these seven gene sets, visualises their overlaps with a SuperVenn diagram, scores each overlap by how unexpected it is, and runs enrichment analysis on the most surprising overlaps.

Worked example. We apply the design to a pair taken from the published DDI–pharmacogenetics knowledge base¹: tacrolimus combined with an azole CYP3A inhibitor (fluconazole or voriconazole), with nephrotoxicity as the adverse event. We pre-specify the expected outputs, the scoring and the decision criteria.

Status of results. This article is a protocol and design paper. No gene sets were downloaded and no overlaps or enrichment statistics were computed for this version. Every statement about what the analysis will find is labelled as a hypothesis. A revision will report the computed results.

Introduction

Polypharmacy is common and rising. A DDI can alter the exposure of the victim drug (a pharmacokinetic, PK, interaction) or change its effect directly (a pharmacodynamic, PD, interaction). Wu et al. trained information-retrieval and information-extraction models on 25 million PubMed abstracts¹. They separated DDI evidence into in-vitro PK, clinical PK and clinical PD studies, and linked DDI pairs to the ADE terms that co-occur with them. Their analysis found 986 DDI pairs supported by all three evidence types and 13,012 pairs supported only by clinical PD evidence. It also generated 119 new CYP3A-related and 18 new CYP2D6-related pharmacogenetic hypotheses¹.

These resources say which pairs interact and which ADEs follow. They rarely say which genes and pathways connect the two drugs to the ADE. Two pieces of public data can fill that gap:

1. Drug perturbation signatures. The LINCS L1000 program profiled transcriptional responses to tens of thousands of compounds across many cell lines. SigCom LINCS serves these as searchable consensus signatures of up- and down-regulated genes².
2. Phenotype-to-gene knowledge. Curated ontologies (the Human Phenotype Ontology, HPO, and the Mammalian Phenotype ontology, MGI-MP), literature co-mention ranking (Geneshot³) and genetic association data (the GWAS Catalog) each give an independent gene set for an adverse event.

If both drugs perturb genes that independently link to the adverse event, those genes are candidate mediators of the interaction. Overlaps that would be unlikely by chance are the most informative. This workflow was designed as a Playbook Workflow Builder (PWB) use case for Common Fund Data Ecosystem (CFDE) resources⁴.

Workflow

The workflow has four stages (Figure 1, given here as a text schematic).

                 DDI pair + adverse event
      (from literature [Wu et al. 2020], FAERS or SIDER)
          /                  |                   \
   Phenotype / ADE         Drug 1               Drug 2
   /      |      \         /    \               /    \
 HPO+   Geneshot  GWAS   L1000  L1000         L1000  L1000
 MGI-MP  lit.     Cat.   down    up           down    up
   \      |      /         \    /               \    /
    \_____|_____/___________\__/_________________\__/
                         |
              SuperVenn of 7 gene sets
                         |
           Score overlaps by "surprise"
                         |
     Enrichment analysis of most surprising overlaps

Figure 1. The "Explain Drug–Drug Interactions (LINCS)" workflow (use case Y1).

Stage 1 — Inputs

• DDI pair and adverse event. One pair (Drug 1, Drug 2) and one ADE term. The sources are the supporting tables of Wu et al.¹ (Table 3 lists substrate–inhibitor pairs with ADEs; Table 4 lists pairs with shared genes), the FDA Adverse Event Reporting System (FAERS) or SIDER.
• Normalisation. Drug names are mapped to the identifiers used in LINCS (for example PubChem CID or LINCS pert_id). ADE terms are mapped from MedDRA preferred terms to the closest HPO and MGI-MP terms. The mapping is recorded for each term.

Stage 2 — Building the seven gene sets

Pre-specified parameters:

• Consensus signatures use the top 250 up- and 250 down-regulated genes.
• Geneshot lists are cut at the top 250 genes by publication count.
• GWAS associations require p < 5 × 10⁻⁸.

All thresholds are reported with the results and varied in a sensitivity analysis (see Stage 3).

Stage 3 — Overlaps and the "surprise" score

The seven sets are drawn as a SuperVenn diagram⁵. Unlike a classic Venn diagram, it shows every non-empty intersection among many sets.

Every gene is classified by (a) which drug-signature sets contain it and (b) how many of the three phenotype sources (sets 5–7) support it. Each drug–phenotype intersection is then scored as follows.

• Over-representation. A one-sided hypergeometric test of the overlap between each drug-signature set and each phenotype set. The background universe is the set of genes measured in the L1000 consensus signatures, not the whole genome. This avoids inflating significance with genes that could never appear in a drug signature.
• Convergence. The number of independent phenotype sources (1–3) that support a gene. A gene supported by HPO/MGI, the literature and GWAS is weighted more heavily than one supported by a single source.
• Interaction specificity. A gene that appears in signatures of both drugs and in at least one phenotype set is a candidate mediator of the interaction rather than of either drug alone. Concordant direction (both drugs down, or both up) suggests an additive effect. Discordant direction suggests opposing effects.
• Null calibration. The observed overlap is compared with 1,000 random drug pairs drawn from LINCS compounds that have the same number of signatures. The empirical percentile is reported alongside the hypergeometric p-value.
• Surprise score. Surprise = −log₁₀(empirical p) × (number of phenotype sources supporting the intersection). Intersections are ranked by this score, and the top intersections pass to Stage 4.
• Sensitivity. Rankings are recomputed with 100 and 500 genes per signature. Only intersections that rank in the top tier at all three cut-offs are called robust.

As the original use case notes, unexpected entries are the point of interest. An overlap containing well-known toxicity genes is reassuring but adds little. An overlap containing genes with no prior link to the ADE but with high scores is a new hypothesis.

Stage 4 — Enrichment of the most surprising overlaps

The genes in the top-ranked intersections go to Enrichr⁶ against pathway (Reactome, WikiPathways, KEGG), Gene Ontology and cell-type libraries. Terms with adjusted p < 0.05 are reported. Terms enriched in the drug–drug–phenotype intersection but not in either single-drug–phenotype intersection are flagged as interaction-specific.

Worked example (pre-specified): tacrolimus + azole antifungal → nephrotoxicity

Why this pair

Wu et al. list tacrolimus as a CYP3A substrate whose interactions with fluconazole and voriconazole are associated with nephrotoxicity, hepatotoxicity and hyperglycaemia. For nephrotoxicity, they cite eleven pharmacogenetic studies supporting a CYP3A link¹. The PK mechanism is established: the azole inhibits CYP3A, so tacrolimus exposure rises. This makes the pair a good positive control. The workflow should recover known calcineurin-inhibitor nephrotoxicity biology, and it can then be asked whether anything beyond raised exposure appears.

A clarithromycin–simvastatin rhabdomyolysis analysis using L1000 data is already on bio.autoxiv. We chose a different pair so that the two analyses complement each other.

Inputs

• Drug 1: tacrolimus (victim, CYP3A substrate).
• Drug 2: fluconazole (perpetrator, CYP3A inhibitor). Voriconazole is used if fluconazole has no adequate L1000 consensus signature.
• ADE: nephrotoxicity (MedDRA). It maps to HPO "Abnormality of the kidney" descendants related to acute kidney injury and tubular injury, and to MGI-MP renal phenotype terms.

Hypotheses (not yet tested)

• H1 — positive control. The tacrolimus signatures will overlap the phenotype sets significantly. The overlap will be enriched for pathways linked to calcineurin-inhibitor nephrotoxicity, such as oxidative stress, TGF-β signalling and fibrosis-related genes.
• H2 — azole contribution. Fluconazole will show smaller overlaps with the nephrotoxicity sets, because its clinical interaction with tacrolimus is mainly pharmacokinetic.
• H3 — interaction-specific genes. If any genes are both (a) in concordant-direction signatures of both drugs and (b) supported by at least two phenotype sources, they are candidate PD mediators. They would point to an interaction component beyond the CYP3A-driven rise in exposure.
• H4 — null outcome. If no drug–drug–phenotype intersection beats the random-pair null, the interaction is best explained by pharmacokinetics alone. That is a useful negative result.

Expected outputs

1. A SuperVenn of the seven gene sets.
2. A ranked table of intersections with hypergeometric p-values, empirical p-values, source counts and surprise scores.
3. Enrichment tables for the top intersections.
4. A gene-level list of candidate mediators with their direction in each drug signature.

Limitations

• Cell-line context. L1000 signatures come mostly from cancer and immortalised cell lines, measured at fixed doses and time points. They may not reflect kidney tissue or clinical exposures. Perturbations that depend on PK, such as the rise in tacrolimus exposure caused by an azole, cannot be seen in single-drug signatures.
• Combinations are not profiled. The workflow uses single-drug signatures as proxies for the combination. True combination effects, whether synergistic or antagonistic, would need dedicated combination profiling.
• Phenotype mapping. Mapping MedDRA terms to HPO and MGI-MP is many-to-many and partly manual. Different mappings can change set 5.
• Literature bias. Geneshot rankings favour well-studied genes. This is partly offset by giving more weight to convergence across sources.
• Source of DDIs. Text-mined DDI–ADE associations are co-occurrence-based and needed manual review in the source study¹. They are inputs to hypothesis generation, not confirmed causal links.
• No results yet. As stated in the Abstract, this version computes nothing. The hypotheses above are untested.

Data and code availability

All inputs are public:

• SigCom LINCS: https://maayanlab.cloud/sigcom-lincs
• Geneshot: https://maayanlab.cloud/geneshot
• GWAS Catalog: https://www.ebi.ac.uk/gwas
• HPO: https://hpo.jax.org
• MGI: https://www.informatics.jax.org
• Enrichr: https://maayanlab.cloud/Enrichr
• Supporting tables from Wu et al.: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7883394/

The workflow can be assembled from PWB building blocks at https://playbook-workflow-builder.cloud/[^4]. Code for the scoring and null calibration will be deposited with the results revision.