An introduction to Refate package
18 December 2025
Package
BiocStyle 2.36.0
Contents
1 Installation
Install the latest development version from GitHub using the devtools package:
library(devtools)
devtools::install() # install the package2 Loading packages and data
suppressMessages(
{
library(Refate)
library(igraph)
library(tidyverse)
library(parallel)
library(networkD3)
}
)
# Load precomputed CPS list and TF-target interaction data
data(CPSs) # Precomputed CPS values for human and mouse
data(TFtarget) # TF-target interaction network
data(drugDB.weight) # Drug-target interaction data
data(TFs_Kinases)
data(StringDB)
data(drugbases_metas)
data(target_data)3 Calculate Cellular Conversion Score (CCS)
# Example log2 fold-change values (replace with your own)
data(fc.fibroblast2hiPSC)
# Calculate CCS for human genes
CCS <- calcCCS(fc.fibroblast2hiPSC, specie = "human")
# View the top CCS scores
head(sort(CCS, decreasing = TRUE))## KRT18 EPCAM SOX2 RBPMS2 ESRP1 KRT19
## 1.5648024 0.9250268 0.7970633 0.6969459 0.6519931 0.6517443# Plot the CCS for known TFs
tf.fibroblast2hiPSC = unique(c("POU5F1", "SOX2", "KLF4", "MYC",
"OCT4", "SOX2", "NANOG", "LIN28", # https://doi.org/10.1126/science.1151526
"ESRRB", "UTF1", "LIN28B", "DPPA2", # https://doi.org/10.1016/j.cell.2012.08.023
"GLIS1", # https://doi.org/10.1038/nature10106
"NR5A2", # https://doi.org/10.1016/j.stem.2009.12.009
"MYCN",
"HMGA1"
))
p <- rank.plot(CCS, gt = tf.fibroblast2hiPSC)
p
4 Prioritize Drugs
# Prioritize drugs
DPS <- calcDrugScore(
CCS,
community_cutoff = 10,
community_cutoff2 = 600,
filtering = TRUE,
drugTarget_cutoff = 500,
TFWeight = 2
)## [1] "The number of sub-graphs is: 43"# View the top-ranked drugs
head(DPS)## ponatinib sorafenib
## 0.9999159 0.9998318
## tinuvin 24morpholinyl8phenyl4h1benzopyran4one
## 0.9997478 0.9996637
## vandetanib wortmannin
## 0.9995796 0.99949555 Build the GRNs
We can build the transcriptional regulatory networks (GRNs) for a given drug. For example, we can build the GRNs for the known drug “Quercetin” by using the following code:
GRNs <- build_GRNs(
drug = "luteolin",
CCS,
kinome,
TFs,
StringDB,
TFtarget,
CCS_threshold = 0.02 # the threshold for selecting genes within the TRN
)
GRNs6 Get drug annotations
We can get drug information from different drug databases
drug_name <- get_drug_name("luteolin")
drug_name ## [1] "Luteolin"drug_annotations <- get_drug_annotations("luteolin")
drug_annotations## $DrugRepurpose.meta
## pert_iname clinical_phase moa target disease_area
## 3524 luteolin Phase 2 glucosidase inhibitor TOP1
## indication direction
## 3524 -1
##
## $CTD.meta
## # A tibble: 207 × 11
## X..ChemicalName ChemicalID CasRN GeneSymbol GeneID GeneForms Organism
## <chr> <chr> <chr> <chr> <int> <chr> <chr>
## 1 Luteolin D047311 491-70-3 ABCB1 5243 protein Homo sapiens
## 2 Luteolin D047311 491-70-3 ABCC1 4363 protein Homo sapiens
## 3 Luteolin D047311 491-70-3 ABCC2 1244 protein Mus musculus
## 4 Luteolin D047311 491-70-3 ABCC4 10257 protein Mus musculus
## 5 Luteolin D047311 491-70-3 ABCG2 9429 protein Homo sapiens
## 6 Luteolin D047311 491-70-3 ACE 1636 protein Oryctolagus …
## 7 Luteolin D047311 491-70-3 ACHE 43 protein Rattus norve…
## 8 Luteolin D047311 491-70-3 ACTA2 59 protein Mus musculus
## 9 Luteolin D047311 491-70-3 AGER 177 mRNA Homo sapiens
## 10 Luteolin D047311 491-70-3 AHR 196 protein Homo sapiens
## # ℹ 197 more rows
## # ℹ 4 more variables: OrganismID <int>, Interaction <chr>,
## # InteractionActions <chr>, PubMedIDs <chr>
##
## $DGIdb.meta
## # A tibble: 32 × 11
## gene_name gene_claim_name entrez_id interaction_claim_sou…¹ interaction_types
## <chr> <chr> <dbl> <chr> <chr>
## 1 ALDH1A1 ALDH1A1 216 DTC <NA>
## 2 ALOX15 ALOX15 246 DTC <NA>
## 3 ALOX15B ALOX15B 247 DTC <NA>
## 4 APEX1 APEX1 328 DTC <NA>
## 5 BAZ2B BAZ2B 29994 DTC <NA>
## 6 BLM BLM 641 DTC <NA>
## 7 CYBB CYBB 1536 DTC <NA>
## 8 CYP1A2 CYP1A2 1544 DTC <NA>
## 9 CYP3A4 CYP3A4 1576 DTC <NA>
## 10 EHMT2 EHMT2 10919 DTC <NA>
## # ℹ 22 more rows
## # ℹ abbreviated name: ¹interaction_claim_source
## # ℹ 6 more variables: drug_claim_name <chr>, drug_claim_primary_name <chr>,
## # drug_name <chr>, drug_concept_id <chr>, interaction_group_score <dbl>,
## # PMIDs <dbl>
##
## $drugBank.meta
## # A tibble: 1 × 15
## primary_key other_keys type created updated name description cas_number
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 DB15584 <NA> small mol… 2019-1… 2020-0… Lute… "" 491-70-3
## # ℹ 7 more variables: unii <chr>, average_mass <chr>, monoisotopic_mass <chr>,
## # state <chr>, synthesis_reference <chr>, fda_label <chr>, msds <chr>drug_description <- get_drug_description("luteolin")
drug_description## [1] ""7 Filtering drugs based on toxicity
filter drugs based on the toxicity information from PubChem database
## we first identify the top drugs
top_drugs <- names(DPS)[which(DPS > 0.9) ]
print(paste0("the number of top drugs: ", length(top_drugs)))## [1] "the number of top drugs: 1189"drug_names <- top_drugs[1:100] # only use the top 100 drugs for demo
data(all_drug_cids)
out <- filter_toxicity(drug_names,
all_drug_cids = all_drug_cids,
json_dir = "tmp/public/PubChem_json")## Parsing PubChem JSON ...## Fetching ChEBI roles ...out$resolution # check which names mapped to which CIDs## # A tibble: 100 × 4
## InputName name_norm CID drug_name
## <chr> <chr> <chr> <chr>
## 1 ponatinib ponatinib 2482… ponatinib
## 2 sorafenib sorafenib 2162… sorafenib
## 3 tinuvin tinuvin 62531 tinuvin
## 4 24morpholinyl8phenyl4h1benzopyran4one 24morpholinyl8phenyl4h… 3973 24morpho…
## 5 vandetanib vandetanib 3081… vandetan…
## 6 wortmannin wortmannin 3121… wortmann…
## 7 dibenzoalpyrene dibenzoalpyrene 9119 dibenzoa…
## 8 asparanina asparanina 2157… asparani…
## 9 melatonin melatonin 896 melatonin
## 10 maneb maneb 3032… maneb
## # ℹ 90 more rowsout$filtered_table %>% dplyr::filter(Keep) %>% dplyr::select(CID, query, MW)## # A tibble: 24 × 3
## CID query MW
## <chr> <chr> <dbl>
## 1 62531 tinuvin 316.
## 2 3973 2-(4-morpholinyl)-8-phenyl-4H-1-benzopyran-4-one 307.
## 3 3973 LY-294002 307.
## 4 896 Melatonin 232.
## 5 176155 SB 203580 377.
## 6 6434217 enzacamene 254.
## 7 5355130 ethylhexyl methoxycinnamate 290.
## 8 5355130 octylmethoxycinnamate 290.
## 9 11167602 regorafenib 483.
## 10 5352624 sulindac sulfide 340.
## # ℹ 14 more rows8 Clutering top drugs based on their target profiles
To characterize mechanistic diversity among top candidate drugs identified from the ESC-to-CNCC conversion analysis, we clustered drugs based on their gene target directionality profiles. This approach groups drugs with similar target engagement patterns and supports the rational selection of compounds spanning distinct mechanisms of action for downstream experimental validation.
data(ESC2CNCC_top_drug)
data(target_data)
data(all_drug_cids)
cids <- all_drug_cids$cid[match(ESC2CNCC_top_drug, all_drug_cids$query)]
target_data_subset <- target_data %>%
dplyr::filter(cid %in% cids) %>%
dplyr::select(cid, genename, direction) %>%
dplyr::filter(!is.na(genename))
mat <- target_data_subset %>%
dplyr::group_by(genename, cid) %>%
dplyr::summarise(direction = direction[1], .groups = "drop") %>% # or mean(direction)
pivot_wider(names_from = cid, values_from = direction, values_fill = 0) %>%
tibble::column_to_rownames("genename") %>%
as.matrix()
colnames(mat) <- ESC2CNCC_top_drug[match(colnames(mat), cids)]
# remove genes with no variation across drugs
mat_filt <- mat[apply(mat, 1, sd) > 0, , drop = FALSE]
mat_filt <- mat_filt[, apply(mat_filt, 2, sd) > 0, drop = FALSE]
# restrict to genes hit by at least k drugs:
k <- 10
mat_filt <- mat_filt[rowSums(mat_filt != 0) >= k, , drop = FALSE]
# Require that each drug hits at least m genes after filtering.
m <- 10
mat_filt2 <- mat_filt[, colSums(mat_filt != 0) >= m, drop = FALSE]
# distance matrix
cor_drug <- cor(mat_filt, method = "pearson", use = "pairwise.complete.obs")
dist_drug <- as.dist(1 - cor_drug)hc <- hclust(dist_drug, method = "average") # average or ward.D2
plot(hc, main = "Drug clustering based on gene targeting", xlab = "", sub = "")
### heatmap visualisation
library(pheatmap)
pheatmap(
t(mat_filt2),
cluster_rows = hc, # <- your precomputed tree
cluster_cols = TRUE, # or another hclust if you have one
scale = "none",
show_rownames = TRUE,
show_colnames = FALSE,
color = colorRampPalette(c("#5074AF", "lightgrey", "#C64032"))(100)
)