TWAS based drug repurposing.

Author

Dr Solal Chauquet

Before you start

This practical is designed to run either offline or on the cluster.
Regardless of your choice, download DrugRepurposing.zip on the CNSG website or on the HPC and unzip it in your directory of choice.

You will end up with the following structure:

DrugRepurposing/
├── DrugRepurposing_prac.html         <- This document
├── data/
│   ├── cmap/                         <- Real CMap signatures (genes x drug-cell-dose)
│   │   ├── CMAP_A375.csv             <- Drug signatures profiled in A375 cell lines
│   │   ├── CMAP_A549.csv             <- Drug signatures profiled in A549 cell lines
|   |   ...
│   │   ├── CMAP_PC3.csv              <- Drug signatures profiled in PC3 cell lines
│   │   └── CMAP_VCAP.csv             <- Drug signatures profiled in VCAP cell lines
│   ├── compound_info.csv             <- Information on the CMap data selected for this practical
│   ├── moa/
│   │   └── MoA_annotations.csv       <- Known mechanism of action for the CMap compounds selected.
│   └── twas/
│   │   ├── sPrediXcan_LDL_Liver.csv  <- real TWAS, LDL cholesterol (GTEx liver eQTL)
│   │   └── sPrediXcan_HDL_Liver.csv  <- real TWAS, HDL cholesterol (GTEx liver eQTL)

Note

For this practical to be feasible within a single hour, both CMap and TWAS data were heavily subsetted. However, all the data used in this practical are real. The CMap signatures used in this practical contain only a few hundred genes and compounds instead of tens of thousands of both. The TWAS files were generated from LDL and HDL¹ and were then subsetted to only contain the most significant genes.

Install the packages we need (only the first time):

install.packages(c("tidyverse", "data.table"))
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("fgsea")

library(tidyverse)
library(data.table)
library(fgsea)

# Set the path to the unzipped directory:
data_dir <- "/path/to/unzip/folder/"
data_dir <- "/Users/uqschauq/Documents/PostDoc/WinterSchool/2026/DrugRepurposing/"

Important

Change the directory in the chunk above to the one where the data was unzipped

Introduction

Bringing a new drug to maket takes well over a decade, costs millions, and only a small fraction of promising drug candidates are ever approved². Drug repurposing allows to side-step this process by finding new indications for already approved, safe drugs. A common often cited example is sildenafil, originally developed for systemic hypertension and now used for erectile dysfunction³, another example includes aspirin, a common pain relief medication now used in both oncology⁴ and cardiovascular prevention⁵.

Previous practicals in this module have explored how disease risk loci identified through GWAS can be mapped to the most likely causal genes, and subsequently to existing drug targets, which can inform repurposing of an approved drug for a new indication. Target-based approaches, such as drug Mendelian Randomisation (drug MR) can be used to validate whether modulating a drug target will have the desired effect on a disease, as well as identify any undesired side effects. While some drugs have a clear target gene/mechanism of action (MoA) (e.g. statins inhibit the HMG-CoA reductase enzyme encoded by the HMGCR gene), for other experimental compounds and approved drugs, their MoA may not be known, and/or they may also have additional unknown targets which we would be unable to test using drug target MR.

When a compound’s MoA/drug target is unknown, another possible approach explored in this practical is to perform a signature-based drug repurposing. Here, we compare gene expression changes associated with a disease of interest to gene expression changes induced by drug exposure. The goal is to identify compounds whose signatures are strongly, negatively correlated with the disease signature, since they exhibit potential to ‘reverse’ the disease gene expression, making them a good candidate for treating the disease of interest.

Signature matching requires only two components:

A disease signature: which genes are up or down-regulated in disease. An in silico approach to generating a disease signature is to conduct a TWAS. Here, GWAS data for the disease of interest is integrated with tissue-specific expression Quantitative Trait Loci (eQTL) data to generate a genetically-predicted disease gene expression signature. Given this is inferred from genetic information, TWAS-identified gene expression association likely points to causal gene expression changes. This offers a major advantage over gene expression measured in case/control studies which can reflect both causal and consequential gene expression changes.
Drug signatures: which genes are up or down-regulated after exposure to a drug. The Connectivity Map (CMap, also called LINCS and L1000)⁶ is a catalog of gene-expression changes induced after exposure to chemical compounds, across a panel of different human cell lines.

In this practical, we will run a TWAS-driven repurposing analysis for two real lipid traits: LDL cholesterol and HDL cholesterol. High LDL-C levels are associated with higher risk of cardiovascular disease. Conversely, higher levels of HDL cholesterol have a protective effect against cardiovascular disease. We will use a TWAS of both LDL and HDL and identify drugs that can lower LDL cholesterol and increase HDL-C.

Transcriptome-wide association study:

Given the liver is the most relevant tissue for lipid metabolism, including cholesterol synthesis, we will conduct a TWAS using a liver-specific eQTL model to impute liver-specific gene expression changes associated with LDL-C and HDL-C.

Let’s load the LDL and HDL TWAS. These are the real S-PrediXcan outputs using liver eQTL models (liver is the principal organ of cholesterol metabolism).

read_twas <- function(path) {
  fread(path) %>%
    # keep genes whose expression model is reliable (Bonferroni on prediction performance)
    filter(pred_perf_pval < 0.05 / n()) %>%
    mutate(metric = zscore) %>%           # we will rank genes by their TWAS z-score
    arrange(zscore) %>%
    select(gene_name, zscore, effect_size, pvalue, metric)
}
ldl <- read_twas(file.path(data_dir, "data/twas", "sPrediXcan_LDL_Liver.csv"))
hdl <- read_twas(file.path(data_dir, "data/twas", "sPrediXcan_HDL_Liver.csv"))
# Look at the most associated genes:
ldl %>% arrange(pvalue) %>% head(10)

The TWAS gives, for every gene, a z-score measuring the association between genetically-predicted expression of that gene and LDL. A positive z-score means that genetically predicted higher expression of the gene is associated with higher LDL.

A negative z-score means that genetically predicted higher expression is associated with lower LDL.

In other words, the sign of the z-score is the direction of effect of that gene’s expression on the trait.

Before performing drug repurposing, we will familiarise ourselves with the TWAS:

Question:

Take the most significant gene in the LDL TWAS. Look it up (GeneCards, OMIM, or just a literature search). What does the protein do? How is it linked to LDL cholesterol biology? Now look at the sign of its z-score — does the direction make sense given what the gene does?
Do the same exercise for the top gene in the HDL TWAS. Are the top genes the same as for LDL, or different? What does that tell you about how specific the genetic signal is to each trait?

Question:

If you wanted a drug to lower LDL, would you look for a positive or a negative correlation between the drug signature and the LDL TWAS? What about a drug to raise HDL?

A drug that reverses the disease signature (downregulating LDL-raising genes and upregulating LDL-lowering genes) will therefore show a negative correlation with the TWAS signature. Such a drug would therefore be a good candidate for repurposing.

Drug signature:

The drug signatures in CMap have been profiled in a panel of different human cell lines. We therefore need to decide which cell line drug signatures we want to query.

You may intuitively think to use the cell line that has the largest number of compounds profiled. Let’s look at the data available:

cell_info <- fread(file.path(data_dir, "/data/compound_info.csv")) %>%
  group_by(cell) %>% summarise(N_available_signature = n())
cell_info %>% arrange(desc(N_available_signature))

Question:

Which cell line has the most compounds profiled? What kind of tissue / disease does that cell line come from? Given that we are studying cholesterol, does that feel like the right tissue?

We will now load the chosen CMap signature matrix (genes in rows, one column per drug-signature) and its metadata.

CELL <- 'CHANGE THIS VALUE TO YOUR CELL'
CELL <- 'MCF7'
fileName <- paste0('CMap_',CELL, '.csv')

cmap     <- fread(file.path(data_dir, "/data/cmap",fileName))
sig_info <- fread(file.path(data_dir, "/data/compound_info.csv")) %>% filter(cell == CELL)
moa      <- fread(file.path(data_dir, "/data/moa/MoA_annotations.csv")) 
ncol(cmap) - 1   # number of drug signatures in this cell line

Important

Select a cell line to use for the repurposing by changing the CELL variable in the chunk above.

Drug-Disease similarity:

For a given TWAS signature, we take the top n_up and bottom n_down genes (ranked by z-score) and compute the Spearman correlation between the TWAS z-scores and each drug’s expression values, over those genes only.

compute_similarity <- function(twas, cmap_matrix, n_up = 25, n_down = 25) {
  twas <- twas %>% arrange(desc(metric))
  sig_genes <- c(head(twas$gene_name, n_up), tail(twas$gene_name, n_down))
  sig <- twas %>% filter(gene_name %in% sig_genes) %>% select(gene_name, metric)

  drug_cols <- setdiff(colnames(cmap_matrix), "gene")
  scores <- vapply(drug_cols, function(d) {
    df <- data.frame(gene_name = cmap_matrix$gene, drug = cmap_matrix[[d]])
    df$drug[is.infinite(df$drug) | is.na(df$drug)] <- 0
    df <- inner_join(sig, df, by = "gene_name")
    suppressWarnings(cor(df$metric, df$drug, method = "spearman"))
  }, numeric(1))
  
  df <- data.frame(colID = drug_cols, score = scores, row.names = NULL)
  df %>% mutate(drug_name = unlist(map(str_split(df$colID, '_'),1)),
                concentration = unlist(map(str_split(df$colID, '_'),4)),
                exposure_duration= unlist(map(str_split(df$colID, '_'),3))) %>%
    select(colID, drug_name, concentration, exposure_duration, score)
}

sim_ldl <- compute_similarity(ldl, cmap, n_up = 25, n_down = 25)
head(sim_ldl[order(sim_ldl$score), ], 5)   # most negative = strongest reversal of the LDL signature

Question:

How many genes were used to calculate the similarity between drug and disease signature?
Given that this number was chosen arbitrarily, how would you choose a number of genes?
Do you think there is any merit in trying multiple gene set sizes with a different number of up and down-regulated genes ?

Let’s look at the most negative correlation for LDL in this cell line:

annotate <- function(sim) {
  sim %>%
    left_join(moa %>% select(-pert_id), by = "drug_name") %>%
    arrange(score)
}

ranked_ldl <- annotate(sim_ldl)
ranked_ldl %>% head(15)

You now have a ranked list of individual compounds. It is tempting to read off the top row and declare it “the candidate”. However, a bit more caution is warranted.

Question:

Most CMap consensus signatures are built from only 3 replicate profiles. If you ran the experiment again, how confident are you that the exact compound sitting in position #1 would still be #1?
Change the number of genes used to calculate the similarity (using the compute_similarity function). Do you have the same top ranked compound? Is the similarity the same?

To explicitly state our point, as most drug signatures are based only on 3 replicates, the robustness of a drug’s gene expression signature is questionable. A compound could be ranked first simply by chance. Robust drug prioritisation requires sensitivity analyses to confirm the top ranking compounds.

Note

While in this practical, we are investigating the similarity of a single gene-set, we would advise using more than one to investigate the stability of the similarity observed, see ⁷ for more information.

Mechanism of action enrichment:

If a particular compound is ranked at the top, we expect that other compounds with the same MoA should also be ranked highly. Therefore, to decrease false discovery and provide greater confidence in our drug prioritisation, we could test whether different compounds that share the same MoA are also ranked highly.

We test this with an enrichment analysis (you can read the GSEA paper to gain a better understanding ⁸). Briefly, we rank all drug signatures by their correlation between drug and TWAS signature, then ask whether the signatures belonging to a given MoA are grouped towards one end rather than scattered at random. If all compounds are in the negative end of the ranked list, then that group is a good target candidate for repurposing.

moa_enrichment <- function(sim, moa_tbl, min_size = 1, nperm = 10000) {
  ranked <- sim %>% filter(!is.na(score)) %>% arrange(score)
  stats <- ranked$score
  names(stats) <- ranked$colID
  stats <- sort(stats, decreasing = TRUE)
  
  
  moa_terms <- ranked %>%
    filter(!is.na(moa), moa != "") %>%
    count(moa) %>% filter(n >= min_size) %>% pull(moa)

  pathways <- lapply(moa_terms, function(m) {
    ranked$colID[ranked$moa == m & !is.na(ranked$moa)]
  })
  names(pathways) <- moa_terms
  pathways <- pathways[lengths(pathways) >= min_size]

  suppressWarnings(
    fgseaSimple(pathways = pathways, stats = stats, nperm = nperm)
  ) %>% as.data.frame() %>% arrange(NES) %>% select(-c(leadingEdge, nMoreExtreme, size))
}

(enr_ldl <- moa_enrichment(ranked_ldl))

Plot every MoA class, ordered by the enrichment score (Normalised enrichment score). A Negative NES means that the class reverses the LDL signature (potential LDL-lowering). Here, we highlight HMGCR inhibitors, a drug class known to decrease LDL cholesterol concentration.

plot_moa <- function(enr, title) {
  enr %>%
    mutate(highlight = ifelse(pathway == "HMGCR inhibitor", "HMGCR inhibitor", "other"),
           pathway = fct_reorder(pathway, NES)) %>%
    ggplot(aes(NES, pathway, fill = highlight)) +
    geom_col() +
    scale_fill_manual(values = c("HMGCR inhibitor" = "darkred", "other" = "grey70")) +
    geom_vline(xintercept = 0, lty = 2) +
    labs(title = title, x = "Normalised Enrichment Score (NES)", y = NULL, fill = NULL) +
    theme_minimal(base_size = 9)
}

plot_moa(enr_ldl, paste0("LDL — MoA enrichment in ", CELL))

# Where do HMGCR inhibitors actually rank?
enr_ldl %>% mutate(rank = row_number()) %>% filter(pathway == "HMGCR inhibitor")

Question:

Look at where HMGCR inhibitors land. Are they among the strongest negative (LDL-reversing) classes? We know statins lower LDL, if they are not at the top here, our analysis has missed a true positive. Why might that be?
Recall the cell line we used and its tissue of origin. CMap signatures are cell-line specific: a drug only produces its cholesterol-relevant transcriptional response in a cell that actually runs the relevant biology. Does it make sense that statins would impact cholesterol linked gene expression in this particular cell line? Form a hypothesis about which cell line would show the statin effect.

Cell line choice:

Cholesterol is synthesised and cleared primarily in the liver, and HMG-CoA reductase (the statin target) is most relevant there. Among the CMap core cell lines, HEPG2 is a hepatocellular (liver) line. It is not the most profiled line, the “more data” is therefore not always the best. When the cell line of choice is unknown, analysis such as the tissue enrichment performed in the FUMA practical can be used to select the cell line of interest.

Let’s switch to the biologically motivated choice and re-run exactly the same analysis.

CELL <- 'HEPG2'
fileName <- paste0('CMap_',CELL, '.csv')
cmap     <- fread(file.path(data_dir, "/data/cmap",fileName))
ncol(cmap) - 1   # number of drug signatures in this cell line

Question:

HEPG2 has fewer signatures than the cell line we used before. In what sense is it nonetheless the better choice? What does this say about choosing analysis parameters by convenience/sample-size versus by biology?

Same code, same parameters (n_up = 25, n_down = 25) — only the cell line has changed.

# calculate the similarity on LDL and HDL:
sim_ldl_hep <- compute_similarity(ldl, cmap, n_up = 25, n_down = 25)
sim_hdl_hep <- compute_similarity(hdl, cmap, n_up = 25, n_down = 25)
# Add the MoA annotation:
ranked_ldl_hep <- annotate(sim_ldl_hep)
ranked_hdl_hep <- annotate(sim_hdl_hep)
# Perform the enrichment calculation:
enr_ldl_hep <- moa_enrichment(ranked_ldl_hep)
enr_hdl_hep <- moa_enrichment(ranked_hdl_hep)
# Plot MoA:
plot_moa(enr_ldl_hep, "LDL — MoA enrichment in HEPG2")
# Get the rank for LDL
enr_ldl_hep %>% mutate(rank = row_number()) %>%
  filter(pathway == "HMGCR inhibitor") %>%
  select(rank, pathway, NES, pval, padj)
# Get the rank for HDL:
enr_hdl_hep %>% mutate(rank = row_number()) %>%
  filter(pathway == "HMGCR inhibitor") %>%
  select(rank, pathway, NES, pval, padj)

Question:

In HEPG2, where do HMGCR inhibitors rank for LDL? Compare the NES and rank to what you saw in the other cell line you used.
Check the sign of the HMGCR NES for LDL versus HDL. Does it match the prediction you made in question 3 (statins reverse LDL → negative; reinforce HDL → positive)?

A repurposing screen is most useful for the hits you didn’t already know about. Now that HEPG2 has enriched HMGCR inhibitors correctly, look at the other MoA classes ranked alongside them.

enr_ldl_hep %>%
  select(pathway, NES, pval, padj) %>% arrange(NES) %>%
  head(10)

Question:

List the MoA classes with the strongest negative (LDL-reversing) enrichment besides HMGCR inhibitors. Pick one and look up its known biology — is there a plausible mechanistic link to cholesterol (e.g. classes acting on lipid uptake, bile-acid handling, or the SREBP pathway)? Would it be worth following up?

Note

This practical is only here as a toy example. While all the data used were real, the CMap matrix has been heavily subsetted. The recovery of the HMGCR effect in biologically relevant cell line only is something we observe in real data, however, the effect of the other mechanisms of action might be due to the subsetting of the data.

Conclusion:

In this practical you ran an end-to-end, genetics-driven drug repurposing analysis on real CMap data, and you established three things:

Starting from nothing but a GWAS-derived TWAS and a database of drug-induced expression signatures, you were able to prioritise HMGCR inhibitors (statins) as a class of drugs for lowering LDL-C. Statins are the first-line treatment for lowering cholesterol and reducing risk of CVD. Statins are also known to moderately increase HDL-C, and this relationship was captured in the weaker positive enrichment of HMGCR inhibitors for HDL-C.
Statins were not prioritised when you used another cell line. The cell line in which drugs are profiled has one of the biggest impacts on drug rankings. If a drug’s target is not expressed in a cell line, the drug signature you observe will be completely different from that measured in a cell line where the target is expressed.
We also saw that individual compound ranks are unreliable when drug signatures are based on a small number of replicates. This is the case for most compounds in CMap. We therefore recommend prioritisation based on drug class enrichment rather than ranking of single compounds when using CMap data. Thinking back to the previous practical, you may also want to run SMR analysis on the target gene based on the MoA of the prioritised drugs for orthogonal validation.

This practical demonstrates how much the choice of drug signature cell line can impact drug prioritisation. However, this is not the only parameter that can impact final drug ranking. The tissue/eQTL model, the similarity metric, the TWAS method and the signature size all influence the final ranking. If you are interested, you can read more about it here⁹

Extension:

You performed the enrichment based on mechanism of action. Is mechanism of action the only possible grouping, or could you think of other group that could be used to perform the enrichment?
This pipeline has prioritised a class of drugs from genetics and transcriptomics. But enrichment is only an observational metric. In the MR practical you learned to estimate the causal effect of a drug target by using genetic variants in/near the target gene as instruments (a “drug-target MR”). Describe how you would follow up your top candidates from today with MR:
- Which gene/target would you instrument for the HMGCR-inhibitor hit?
- What exposure and outcome would you use?
- What would a significant MR result add that the CMap enrichment alone cannot?
- How could MR also flag potential on-target side effects of a candidate before any wet-lab work?

References:

Footnotes

Graham, S. E. et al. The power of genetic diversity in genome-wide association studies of lipids. Nature 600, 675–679 (2021).↩︎
Yamaguchi, S., Kaneko, M. & Narukawa, M. Approval success rates of drug candidates based on target, action, modality, application, and their combinations. Clin. Transl. Sci. 14, 1113–1122 (2021).↩︎
Cruz-Burgos, M. et al. New Approaches in Oncology for Repositioning Drugs: The Case of PDE5 Inhibitor Sildenafil. Front. Oncol. 11, 627229 (2021).↩︎
Sun, M. et al. Role of aspirin in cancer prevention. Cancer Treat. Res. Commun. 43, 100884 (2025).↩︎
Cofer, L. B., Barrett, T. J. & Berger, J. S. Aspirin for the Primary Prevention of Cardiovascular Disease: Time for a Platelet-Guided Approach. Arterioscler. Thromb. Vasc. Biol. 42, 1207–1216 (2022).↩︎
Subramanian, A. et al. A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell 171, 1437–1452.e17 (2017).↩︎
Chauquet, S. et al. (2026). [Drug repurposing from transcriptome-wide association studies]. medRxiv 2026.04.22.26349470.↩︎
Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550 (2005).↩︎
Chauquet, S. et al. (2026). [Drug repurposing from transcriptome-wide association studies]. medRxiv 2026.04.22.26349470.↩︎

--- title: "TWAS based drug repurposing." author: "Dr Solal Chauquet" format: html: toc: true toc-location: left code-fold: false code-tools: true theme: cosmo embed-resources: true editor: visual --- ## Before you start This practical is designed to run either offline or on the cluster.\ Regardless of your choice, download `DrugRepurposing.zip` on the CNSG website or on the HPC and unzip it in your directory of choice.\ You will end up with the following structure: ``` DrugRepurposing/ ├── DrugRepurposing_prac.html <- This document ├── data/ │ ├── cmap/ <- Real CMap signatures (genes x drug-cell-dose) │ │ ├── CMAP_A375.csv <- Drug signatures profiled in A375 cell lines │ │ ├── CMAP_A549.csv <- Drug signatures profiled in A549 cell lines | | ... │ │ ├── CMAP_PC3.csv <- Drug signatures profiled in PC3 cell lines │ │ └── CMAP_VCAP.csv <- Drug signatures profiled in VCAP cell lines │ ├── compound_info.csv <- Information on the CMap data selected for this practical │ ├── moa/ │ │ └── MoA_annotations.csv <- Known mechanism of action for the CMap compounds selected. │ └── twas/ │ │ ├── sPrediXcan_LDL_Liver.csv <- real TWAS, LDL cholesterol (GTEx liver eQTL) │ │ └── sPrediXcan_HDL_Liver.csv <- real TWAS, HDL cholesterol (GTEx liver eQTL) ``` ::: callout-note For this practical to be feasible within a single hour, both CMap and TWAS data were heavily subsetted. However, all the data used in this practical are *real*. The CMap signatures used in this practical contain only a few hundred genes and compounds instead of tens of thousands of both. The TWAS files were generated from LDL and HDL[^1] and were then subsetted to only contain the most significant genes. ::: [^1]: Graham, S. E. et al. The power of genetic diversity in genome-wide association studies of lipids. Nature 600, 675–679 (2021). Install the packages we need (only the first time): ```{r} #| eval: false install.packages(c("tidyverse", "data.table")) if (!require("BiocManager")) install.packages("BiocManager") BiocManager::install("fgsea") ``` ```{r} #| eval: false #| message: false #| warning: false library(tidyverse) library(data.table) library(fgsea) # Set the path to the unzipped directory: data_dir <- "/path/to/unzip/folder/" data_dir <- "/Users/uqschauq/Documents/PostDoc/WinterSchool/2026/DrugRepurposing/" ``` ::: callout-important Change the directory in the chunk above to the one where the data was unzipped ::: ## Introduction Bringing a new drug to maket takes well over a decade, costs millions, and only a small fraction of promising drug candidates are ever approved[^2]. Drug repurposing allows to side-step this process by finding new indications for already approved, safe drugs. A common often cited example is sildenafil, originally developed for systemic hypertension and now used for erectile dysfunction[^3], another example includes aspirin, a common pain relief medication now used in both oncology[^4] and cardiovascular prevention[^5]. [^2]: Yamaguchi, S., Kaneko, M. & Narukawa, M. Approval success rates of drug candidates based on target, action, modality, application, and their combinations. *Clin. Transl. Sci.* **14**, 1113–1122 (2021). [^3]: Cruz-Burgos, M. *et al.* New Approaches in Oncology for Repositioning Drugs: The Case of PDE5 Inhibitor Sildenafil. *Front. Oncol.* **11**, 627229 (2021). [^4]: Sun, M. et al. Role of aspirin in cancer prevention. Cancer Treat. Res. Commun. 43, 100884 (2025). [^5]: Cofer, L. B., Barrett, T. J. & Berger, J. S. Aspirin for the Primary Prevention of Cardiovascular Disease: Time for a Platelet-Guided Approach. Arterioscler. Thromb. Vasc. Biol. 42, 1207–1216 (2022). Previous practicals in this module have explored how disease risk loci identified through GWAS can be mapped to the most likely causal genes, and subsequently to existing drug targets, which can inform repurposing of an approved drug for a new indication. Target-based approaches, such as drug Mendelian Randomisation (drug MR) can be used to validate whether modulating a drug target will have the desired effect on a disease, as well as identify any undesired side effects. While some drugs have a clear target gene/mechanism of action (MoA) (e.g. statins inhibit the HMG-CoA reductase enzyme encoded by the HMGCR gene), for other experimental compounds and approved drugs, their MoA may not be known, and/or they may also have additional unknown targets which we would be unable to test using drug target MR. When a compound's MoA/drug target is unknown, another possible approach explored in this practical is to perform a signature-based drug repurposing. Here, we compare gene expression changes associated with a disease of interest to gene expression changes induced by drug exposure. The goal is to identify compounds whose signatures are strongly, negatively correlated with the disease signature, since they exhibit potential to 'reverse' the disease gene expression, making them a good candidate for treating the disease of interest. Signature matching requires only two components: 1. **A disease signature**: which genes are up or down-regulated in disease. An in silico approach to generating a disease signature is to conduct a TWAS. Here, GWAS data for the disease of interest is integrated with tissue-specific expression Quantitative Trait Loci (eQTL) data to generate a genetically-predicted disease gene expression signature. Given this is inferred from genetic information, TWAS-identified gene expression association likely points to causal gene expression changes. This offers a major advantage over gene expression measured in case/control studies which can reflect both causal and consequential gene expression changes. 2. **Drug signatures**: which genes are up or down-regulated after exposure to a drug. The **Connectivity Map** (CMap, also called LINCS and L1000)[^6] is a catalog of gene-expression changes induced after exposure to chemical compounds, across a panel of different human cell lines. [^6]: Subramanian, A. *et al.* A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. *Cell* **171**, 1437–1452.e17 (2017). In this practical, we will run a TWAS-driven repurposing analysis for two real lipid traits: **LDL cholesterol** and **HDL cholesterol**. High LDL-C levels are associated with higher risk of cardiovascular disease. Conversely, higher levels of HDL cholesterol have a protective effect against cardiovascular disease. We will use a TWAS of both LDL and HDL and identify drugs that can lower LDL cholesterol and increase HDL-C. ## Transcriptome-wide association study: Given the liver is the most relevant tissue for lipid metabolism, including cholesterol synthesis, we will conduct a TWAS using a liver-specific eQTL model to impute liver-specific gene expression changes associated with LDL-C and HDL-C. Let's load the LDL and HDL TWAS. These are the **real** S-PrediXcan outputs using **liver** eQTL models (liver is the principal organ of cholesterol metabolism). ```{r} #| eval: false read_twas <- function(path) { fread(path) %>% # keep genes whose expression model is reliable (Bonferroni on prediction performance) filter(pred_perf_pval < 0.05 / n()) %>% mutate(metric = zscore) %>% # we will rank genes by their TWAS z-score arrange(zscore) %>% select(gene_name, zscore, effect_size, pvalue, metric) } ldl <- read_twas(file.path(data_dir, "data/twas", "sPrediXcan_LDL_Liver.csv")) hdl <- read_twas(file.path(data_dir, "data/twas", "sPrediXcan_HDL_Liver.csv")) # Look at the most associated genes: ldl %>% arrange(pvalue) %>% head(10) ``` The TWAS gives, for every gene, a **z-score** measuring the association between genetically-predicted expression of that gene and LDL. A **positive** z-score means that genetically predicted *higher* expression of the gene is associated with *higher* LDL. A **negative** z-score means that genetically predicted *higher* expression is associated with *lower* LDL. In other words, the sign of the z-score is the **direction of effect** of that gene's expression on the trait. Before performing drug repurposing, we will familiarise ourselves with the TWAS: ::: callout-tip ## Question: 1. Take the most significant gene in the LDL TWAS. Look it up (GeneCards, OMIM, or just a literature search). *What does the protein do? How is it linked to LDL cholesterol biology?* Now look at the **sign** of its z-score — does the direction make sense given what the gene does? 2. Do the same exercise for the top gene in the **HDL** TWAS. Are the top genes the same as for LDL, or different? What does that tell you about how specific the genetic signal is to each trait? ::: ::: callout-tip ## Question: 3. If you wanted a drug to *lower* LDL, would you look for a *positive* or a *negative* correlation between the drug signature and the LDL TWAS? What about a drug to *raise* HDL? ::: A drug that **reverses** the disease signature (downregulating LDL-raising genes and upregulating LDL-lowering genes) will therefore show a negative correlation with the TWAS signature. Such a drug would therefore be a good candidate for repurposing. ## Drug signature: The drug signatures in CMap have been profiled in a panel of different human cell lines. We therefore need to decide which cell line drug signatures we want to query. You may intuitively think to use the cell line that has the largest number of compounds profiled. Let’s look at the data available: ```{r} #| eval: false cell_info <- fread(file.path(data_dir, "/data/compound_info.csv")) %>% group_by(cell) %>% summarise(N_available_signature = n()) cell_info %>% arrange(desc(N_available_signature)) ``` ::: callout-tip ## Question: 4. Which cell line has the most compounds profiled? What kind of tissue / disease does that cell line come from? Given that we are studying **cholesterol**, does that feel like the right tissue? ::: We will now load the chosen CMap signature matrix (genes in rows, one column per drug-signature) and its metadata. ```{r} #| eval: false CELL <- 'CHANGE THIS VALUE TO YOUR CELL' CELL <- 'MCF7' fileName <- paste0('CMap_',CELL, '.csv') cmap <- fread(file.path(data_dir, "/data/cmap",fileName)) sig_info <- fread(file.path(data_dir, "/data/compound_info.csv")) %>% filter(cell == CELL) moa <- fread(file.path(data_dir, "/data/moa/MoA_annotations.csv")) ncol(cmap) - 1 # number of drug signatures in this cell line ``` ::: callout-important Select a cell line to use for the repurposing by changing the CELL variable in the chunk above. ::: ## Drug-Disease similarity: For a given TWAS signature, we take the top `n_up` and bottom `n_down` genes (ranked by z-score) and compute the **Spearman correlation** between the TWAS z-scores and each drug's expression values, over those genes only. ```{r} #| eval: false compute_similarity <- function(twas, cmap_matrix, n_up = 25, n_down = 25) { twas <- twas %>% arrange(desc(metric)) sig_genes <- c(head(twas$gene_name, n_up), tail(twas$gene_name, n_down)) sig <- twas %>% filter(gene_name %in% sig_genes) %>% select(gene_name, metric) drug_cols <- setdiff(colnames(cmap_matrix), "gene") scores <- vapply(drug_cols, function(d) { df <- data.frame(gene_name = cmap_matrix$gene, drug = cmap_matrix[[d]]) df$drug[is.infinite(df$drug) | is.na(df$drug)] <- 0 df <- inner_join(sig, df, by = "gene_name") suppressWarnings(cor(df$metric, df$drug, method = "spearman")) }, numeric(1)) df <- data.frame(colID = drug_cols, score = scores, row.names = NULL) df %>% mutate(drug_name = unlist(map(str_split(df$colID, '_'),1)), concentration = unlist(map(str_split(df$colID, '_'),4)), exposure_duration= unlist(map(str_split(df$colID, '_'),3))) %>% select(colID, drug_name, concentration, exposure_duration, score) } sim_ldl <- compute_similarity(ldl, cmap, n_up = 25, n_down = 25) head(sim_ldl[order(sim_ldl$score), ], 5) # most negative = strongest reversal of the LDL signature ``` ::: callout-tip ## Question: 5. How many genes were used to calculate the similarity between drug and disease signature? 6. Given that this number was chosen arbitrarily, how would you choose a number of genes? 7. Do you think there is any merit in trying multiple gene set sizes with a different number of up and down-regulated genes ? ::: Let's look at the most negative correlation for LDL in this cell line: ```{r} #| eval: false annotate <- function(sim) { sim %>% left_join(moa %>% select(-pert_id), by = "drug_name") %>% arrange(score) } ranked_ldl <- annotate(sim_ldl) ranked_ldl %>% head(15) ``` You now have a ranked list of individual compounds. It is tempting to read off the top row and declare it "the candidate". However, a bit more caution is warranted. ::: callout-tip ## Question: 8. Most CMap consensus signatures are built from only **3** replicate profiles. If you ran the experiment again, how confident are you that the exact compound sitting in position #1 would still be #1? 9. Change the number of genes used to calculate the similarity (using the `compute_similarity` function). Do you have the same top ranked compound? Is the similarity the same? ::: To explicitly state our point, as most drug signatures are based only on 3 replicates, the robustness of a drug's gene expression signature is questionable. A compound could be ranked first simply by chance. Robust drug prioritisation requires sensitivity analyses to confirm the top ranking compounds. ::: callout-note While in this practical, we are investigating the similarity of a single gene-set, we would advise using more than one to investigate the stability of the similarity observed, see [^7] for more information. ::: [^7]: Chauquet, S. *et al.* (2026). \[Drug repurposing from transcriptome-wide association studies\]. *medRxiv* 2026.04.22.26349470. ## Mechanism of action enrichment: If a particular compound is ranked at the top, we expect that other compounds with the same MoA should also be ranked highly. Therefore, to decrease false discovery and provide greater confidence in our drug prioritisation, we could test whether different compounds that share the same MoA are also ranked highly. We test this with an **enrichment** analysis (you can read the GSEA paper to gain a better understanding [^8]). Briefly, we rank all drug signatures by their correlation between drug and TWAS signature, then ask whether the signatures belonging to a given MoA are grouped towards one end rather than scattered at random. If all compounds are in the negative end of the ranked list, then that group is a good target candidate for repurposing. [^8]: Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550 (2005). ```{r} #| eval: false moa_enrichment <- function(sim, moa_tbl, min_size = 1, nperm = 10000) { ranked <- sim %>% filter(!is.na(score)) %>% arrange(score) stats <- ranked$score names(stats) <- ranked$colID stats <- sort(stats, decreasing = TRUE) moa_terms <- ranked %>% filter(!is.na(moa), moa != "") %>% count(moa) %>% filter(n >= min_size) %>% pull(moa) pathways <- lapply(moa_terms, function(m) { ranked$colID[ranked$moa == m & !is.na(ranked$moa)] }) names(pathways) <- moa_terms pathways <- pathways[lengths(pathways) >= min_size] suppressWarnings( fgseaSimple(pathways = pathways, stats = stats, nperm = nperm) ) %>% as.data.frame() %>% arrange(NES) %>% select(-c(leadingEdge, nMoreExtreme, size)) } (enr_ldl <- moa_enrichment(ranked_ldl)) ``` Plot every MoA class, ordered by the enrichment score (Normalised enrichment score). A Negative NES means that the class reverses the LDL signature (potential LDL-lowering). Here, we highlight HMGCR inhibitors, a drug class known to decrease LDL cholesterol concentration. ```{r} #| eval: false #| fig-width: 7 #| fig-height: 9 plot_moa <- function(enr, title) { enr %>% mutate(highlight = ifelse(pathway == "HMGCR inhibitor", "HMGCR inhibitor", "other"), pathway = fct_reorder(pathway, NES)) %>% ggplot(aes(NES, pathway, fill = highlight)) + geom_col() + scale_fill_manual(values = c("HMGCR inhibitor" = "darkred", "other" = "grey70")) + geom_vline(xintercept = 0, lty = 2) + labs(title = title, x = "Normalised Enrichment Score (NES)", y = NULL, fill = NULL) + theme_minimal(base_size = 9) } plot_moa(enr_ldl, paste0("LDL — MoA enrichment in ", CELL)) # Where do HMGCR inhibitors actually rank? enr_ldl %>% mutate(rank = row_number()) %>% filter(pathway == "HMGCR inhibitor") ``` ::: callout-tip ## Question: 10. Look at where **HMGCR inhibitors** land. Are they among the strongest *negative* (LDL-reversing) classes? We *know* statins lower LDL, if they are **not** at the top here, our analysis has missed a true positive. Why might that be? 11. Recall the cell line we used and its tissue of origin. CMap signatures are **cell-line specific**: a drug only produces its cholesterol-relevant transcriptional response in a cell that actually runs the relevant biology. Does it make sense that statins would impact cholesterol linked gene expression in *this* particular cell line? Form a hypothesis about which cell line *would* show the statin effect. ::: ## Cell line choice: Cholesterol is synthesised and cleared primarily in the **liver**, and HMG-CoA reductase (the statin target) is most relevant there. Among the CMap core cell lines, **HEPG2** is a hepatocellular (liver) line. It is *not* the most profiled line, the "more data" is therefore not always the best. When the cell line of choice is unknown, analysis such as the tissue enrichment performed in the FUMA practical can be used to select the cell line of interest. Let's switch to the biologically motivated choice and re-run **exactly the same analysis**. ```{r} #| eval: false CELL <- 'HEPG2' fileName <- paste0('CMap_',CELL, '.csv') cmap <- fread(file.path(data_dir, "/data/cmap",fileName)) ncol(cmap) - 1 # number of drug signatures in this cell line ``` ::: callout-tip ## Question: 12. HEPG2 has *fewer* signatures than the cell line we used before. In what sense is it nonetheless the *better* choice? What does this say about choosing analysis parameters by convenience/sample-size versus by biology? ::: Same code, same parameters (`n_up = 25`, `n_down = 25`) — only the cell line has changed. ```{r} #| eval: false # calculate the similarity on LDL and HDL: sim_ldl_hep <- compute_similarity(ldl, cmap, n_up = 25, n_down = 25) sim_hdl_hep <- compute_similarity(hdl, cmap, n_up = 25, n_down = 25) # Add the MoA annotation: ranked_ldl_hep <- annotate(sim_ldl_hep) ranked_hdl_hep <- annotate(sim_hdl_hep) # Perform the enrichment calculation: enr_ldl_hep <- moa_enrichment(ranked_ldl_hep) enr_hdl_hep <- moa_enrichment(ranked_hdl_hep) # Plot MoA: plot_moa(enr_ldl_hep, "LDL — MoA enrichment in HEPG2") # Get the rank for LDL enr_ldl_hep %>% mutate(rank = row_number()) %>% filter(pathway == "HMGCR inhibitor") %>% select(rank, pathway, NES, pval, padj) # Get the rank for HDL: enr_hdl_hep %>% mutate(rank = row_number()) %>% filter(pathway == "HMGCR inhibitor") %>% select(rank, pathway, NES, pval, padj) ``` ::: callout-tip ## Question: 13. In HEPG2, where do HMGCR inhibitors rank for **LDL**? Compare the NES and rank to what you saw in the other cell line you used. 14. Check the **sign** of the HMGCR NES for LDL versus HDL. Does it match the prediction you made in question 3 (statins *reverse* LDL → negative; *reinforce* HDL → positive)? ::: A repurposing screen is most useful for the hits you *didn't* already know about. Now that HEPG2 has enriched HMGCR inhibitors correctly, look at the **other MoA classes ranked alongside them**. ```{r} #| eval: false enr_ldl_hep %>% select(pathway, NES, pval, padj) %>% arrange(NES) %>% head(10) ``` ::: callout-tip ## Question: 15. List the MoA classes with the strongest negative (LDL-reversing) enrichment besides HMGCR inhibitors. Pick one and look up its known biology — is there a plausible mechanistic link to cholesterol (e.g. classes acting on lipid uptake, bile-acid handling, or the SREBP pathway)? Would it be worth following up? ::: ::: callout-note This practical is only here as a toy example. While all the data used were real, the CMap matrix has been heavily subsetted. The recovery of the HMGCR effect in biologically relevant cell line only is something we observe in real data, however, the effect of the other mechanisms of action might be due to the subsetting of the data. ::: ## Conclusion: In this practical you ran an end-to-end, genetics-driven drug repurposing analysis on real CMap data, and you established three things: 1. Starting from nothing but a GWAS-derived TWAS and a database of drug-induced expression signatures, you were able to prioritise HMGCR inhibitors (statins) as a class of drugs for lowering LDL-C. Statins are the first-line treatment for lowering cholesterol and reducing risk of CVD. Statins are also known to moderately increase HDL-C, and this relationship was captured in the weaker positive enrichment of HMGCR inhibitors for HDL-C. 2. Statins were not prioritised when you used another cell line. The cell line in which drugs are profiled has one of the biggest impacts on drug rankings. If a drug's target is not expressed in a cell line, the drug signature you observe will be completely different from that measured in a cell line where the target is expressed. 3. We also saw that individual compound ranks are unreliable when drug signatures are based on a small number of replicates. This is the case for most compounds in CMap. We therefore recommend prioritisation based on drug class enrichment rather than ranking of single compounds when using CMap data. Thinking back to the previous practical, you may also want to run SMR analysis on the target gene based on the MoA of the prioritised drugs for orthogonal validation. This practical demonstrates how much the choice of drug signature cell line can impact drug prioritisation. However, this is not the only parameter that can impact final drug ranking. The tissue/eQTL model, the similarity metric, the TWAS method and the signature size all influence the final ranking. If you are interested, you can read more about it here[^9] [^9]: Chauquet, S. *et al.* (2026). \[Drug repurposing from transcriptome-wide association studies\]. *medRxiv* 2026.04.22.26349470. ## Extension: ::: callout-tip ## Extension: - You performed the enrichment based on mechanism of action. Is mechanism of action the only possible grouping, or could you think of other group that could be used to perform the enrichment? - This pipeline has *prioritised* a class of drugs from genetics and transcriptomics. But enrichment is only an observational metric. In the MR practical you learned to estimate the causal effect of a drug target by using genetic variants in/near the target gene as instruments (a "drug-target MR"). Describe how you would follow up your top candidates from today with MR: - Which gene/target would you instrument for the HMGCR-inhibitor hit? - What exposure and outcome would you use? - What would a significant MR result add that the CMap enrichment alone cannot? - How could MR also flag potential **on-target side effects** of a candidate before any wet-lab work? ::: ## References: