Skip to content

Add MPNST treated drug-combination experiments #488

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jjacobson95 wants to merge 1 commit into
base: mpnst-treated-samples-and-omics
Choose a base branch
from
Open

Add MPNST treated drug-combination experiments #488

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
284 changes: 235 additions & 49 deletions coderbuild/mpnst/02_get_drug_data.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,48 @@
# 02_get_drug_data.R
#!/usr/bin/env Rscript

# Combined Drug List Extraction for MPNST & MPNST‑PDX
# Combined Drug List Extraction for MPNST & MPNST-PDX

library(data.table)
library(dplyr)
library(stringr)
library(synapser)
library(reticulate)

# ----------------------------
# DEBUG helper
# ----------------------------
debug_file_summary <- function(path, label) {
message("\n========== DEBUG: ", label, " ==========")
message("Path: ", path)
if (is.null(path) || is.na(path) || !file.exists(path)) {
message("File missing.")
return(invisible(NULL))
}

df <- tryCatch(fread(path, sep="\t", header=TRUE), error=function(e) {
message("fread failed: ", conditionMessage(e))
return(NULL)
})
if (is.null(df)) return(invisible(NULL))

message("Rows: ", nrow(df), " Cols: ", ncol(df))
message("Columns: ", paste(names(df), collapse=", "))

if ("chem_name" %in% names(df)) {
cn <- unique(tolower(trimws(as.character(df$chem_name))))
message("Unique chem_name (first 25): ", paste(head(cn, 25), collapse=", "))
message("Has chem_name == 'verteporfin'? ", any(cn == "verteporfin", na.rm = TRUE))
message("Has chem_name == 'wu713d62n9'? ", any(cn == "wu713d62n9", na.rm = TRUE))
}
if ("pubchem_id" %in% names(df)) {
pc <- unique(trimws(as.character(df$pubchem_id)))
message("Unique pubchem_id (first 10): ", paste(head(pc, 10), collapse=", "))
}

invisible(df)
}

# 0) Args & login
args <- commandArgs(trailingOnly = TRUE)
if (length(args) < 1) {
Expand All @@ -25,38 +60,54 @@ synLogin(authToken = token)
manifest <- synTableQuery("select * from syn53503360")$asDataFrame() %>%
rename(common_name = Sample)

# 2) PDX‑sourced drugs via annotations
pdx_df <- manifest %>%
select(common_name, PDX_Drug_Data) %>%
distinct() %>%
filter(!is.na(PDX_Drug_Data))
# 2) PDX-sourced drugs via Synapse IDs stored in the manifest drug column.
# Column was renamed PDX_Drug_Data → PDXDrugData in 2025; values may now be
# file names rather than Synapse IDs, so filter to syn* entries only.
pdx_col <- if ("PDXDrugData" %in% names(manifest)) "PDXDrugData" else
if ("PDX_Drug_Data" %in% names(manifest)) "PDX_Drug_Data" else NULL
pdx_drugs <- character(0)
if (!is.null(pdx_col)) {
pdx_df <- manifest %>%
select(common_name, all_of(pdx_col)) %>%
rename(pdx_data = all_of(pdx_col)) %>%
distinct() %>%
filter(!is.na(pdx_data))

pdx_ids <- unique(unlist(strsplit(pdx_df$PDX_Drug_Data, ",")))
pdx_ids <- pdx_ids[ pdx_ids != "" & !is.na(pdx_ids) & pdx_ids != "NA" ]

get_pdx_drugs <- function(synid) {
# Query the metadata table for this file's experimentalCondition
q <- sprintf(
"select experimentalCondition from syn21993642 where id='%s'",
synid
)
df <- synTableQuery(q)$asDataFrame()
if (nrow(df)==0) return(character(0))
# Split on semicolon, lowercase and drop empties
conds <- unlist(strsplit(df$experimentalCondition, ";"))
tolower(conds[conds!=""])
}
pdx_ids <- unique(unlist(lapply(pdx_df$pdx_data, function(x) {
strsplit(paste(trimws(unlist(x)), collapse=","), ",")[[1]]
})))
pdx_ids <- pdx_ids[ pdx_ids != "" & !is.na(pdx_ids) & pdx_ids != "NA" ]
pdx_ids <- pdx_ids[ grepl("^syn[0-9]+$", trimws(pdx_ids), ignore.case = TRUE) ]

pdx_drugs <- unique(unlist(lapply(pdx_ids, get_pdx_drugs)))
pdx_drugs <- setdiff(pdx_drugs, "control")
get_pdx_drugs <- function(synid) {
q <- sprintf(
"select experimentalCondition from syn21993642 where id='%s'",
synid
)
df <- synTableQuery(q)$asDataFrame()
if (nrow(df)==0) return(character(0))
conds <- unlist(strsplit(as.character(df$experimentalCondition), ";"))
tolower(conds[conds!=""])
}

if (length(pdx_ids) > 0) {
pdx_drugs <- unique(unlist(lapply(pdx_ids, get_pdx_drugs)))
pdx_drugs <- setdiff(pdx_drugs, "control")
} else {
message("No Synapse IDs found in ", pdx_col, "; skipping PDX Synapse drug extraction.")
}
} else {
message("No PDX drug column found in manifest; skipping PDX drug extraction.")
}

# 3) MicroTissue‑sourced drugs via table "children"
# 3) MicroTissue-sourced drugs via table "children"
mts_df <- manifest %>%
select(common_name, MicroTissueDrugFolder) %>%
filter(!is.na(MicroTissueDrugFolder))

mts_ids <- unique(unlist(strsplit(mts_df$MicroTissueDrugFolder, ",")))
mts_ids <- unique(unlist(lapply(mts_df$MicroTissueDrugFolder, function(x) {
strsplit(paste(trimws(unlist(x)), collapse=","), ",")[[1]]
})))
mts_ids <- mts_ids[mts_ids != "" & !is.na(mts_ids) & mts_ids != "NA"]

get_mts_drugs <- function(parentId) {
Expand All @@ -69,17 +120,93 @@ get_mts_drugs <- function(parentId) {

mts_drugs <- unique(unlist(lapply(mts_ids, get_mts_drugs)))

# 3.5) IMPROVE-sourced drugs via improve_sample_id (syn69801348)
improve_sample_drugs <- character(0)
improve_path <- tryCatch(synGet("syn69801348")$path, error=function(e) NA)
if (!is.na(improve_path) && file.exists(improve_path)) {
improve_tab <- fread(improve_path, sep="\t", header=TRUE, fill=TRUE)
if ("improve_sample_id" %in% names(improve_tab)) {
ids <- improve_tab$improve_sample_id %>% as.character() %>% unique()
ids <- ids[!is.na(ids) & ids != ""]
drug_part <- sub("^.*_", "", ids)
drug_tokens <- unlist(strsplit(drug_part, "\\+"))
improve_sample_drugs <- unique(tolower(trimws(drug_tokens)))
improve_sample_drugs <- improve_sample_drugs[improve_sample_drugs != ""]
}
}

# 3.6) Add drugs from treatment lookup list (syn64608501) Full_name column
lookup_drugs <- character(0)
lookup_path <- tryCatch(synGet("syn64608501")$path, error=function(e) NA)
if (!is.na(lookup_path) && file.exists(lookup_path)) {
lookup_dt <- tryCatch(
fread(lookup_path, sep = "\t", header = TRUE, fill = TRUE),
error = function(e) NULL
)

if (is.null(lookup_dt)) {
lookup_dt <- tryCatch(
as.data.table(readxl::read_excel(lookup_path)),
error = function(e) NULL
)
}

if (!is.null(lookup_dt)) {
full_col <- NULL
if ("Full_name" %in% names(lookup_dt)) full_col <- "Full_name"
if (is.null(full_col) && "Full Name" %in% names(lookup_dt)) full_col <- "Full Name"

if (!is.null(full_col)) {
lookup_drugs <- unique(tolower(trimws(as.character(lookup_dt[[full_col]]))))
lookup_drugs <- lookup_drugs[!is.na(lookup_drugs) & nzchar(lookup_drugs)]
lookup_drugs <- setdiff(lookup_drugs, c("control", "dimethyl_sulfoxide", "dimethyl sulfoxide", "dmso"))
} else {
message("\n========== DEBUG: DRUG LOOKUP LIST (syn64608501) ==========")
message("WARNING: Could not find Full_name column in lookup file. Columns: ",
paste(names(lookup_dt), collapse = ", "))
}
} else {
message("\n========== DEBUG: DRUG LOOKUP LIST (syn64608501) ==========")
message("WARNING: Could not read lookup file at: ", lookup_path)
}
} else {
message("\n========== DEBUG: DRUG LOOKUP LIST (syn64608501) ==========")
message("WARNING: lookup file missing/unavailable for syn64608501")
}

message("\n========== DEBUG: DRUG LOOKUP LIST (syn64608501) ==========")
message("lookup_path: ", ifelse(is.na(lookup_path), "<NA>", lookup_path))
message("lookup_drugs count: ", length(lookup_drugs))
if (length(lookup_drugs) > 0) {
message("First 40 lookup_drugs: ", paste(head(lookup_drugs, 40), collapse = ", "))
}

# 4) Combine and fix bad names
all_drugs <- unique(c(pdx_drugs, mts_drugs))
all_drugs <- unique(c(pdx_drugs, mts_drugs, improve_sample_drugs, lookup_drugs))
all_drugs[all_drugs == "pd901"] <- "pd-0325901"
# message("Combined drug list: ", paste(all_drugs, collapse=", "))
all_drugs[all_drugs == "rmc4630"] <- "rmc-4630"

# IMPORTANT: drop NA/empty
all_drugs <- all_drugs[!is.na(all_drugs) & nzchar(all_drugs)]

# ---- NEW: replace verteporfin with WU713D62N9 so we can do ONE name-based call
# (keep exact casing you provided)
if (any(all_drugs == "verteporfin", na.rm = TRUE)) {
all_drugs[all_drugs == "verteporfin"] <- "WU713D62N9"
}

# 5) Read old‑drug files or initialize empty
message("\n========== DEBUG: DRUG LISTS ==========")
message("all_drugs count: ", length(all_drugs))
message("Contains 'verteporfin' in all_drugs? ", any(all_drugs == "verteporfin", na.rm=TRUE))
message("Contains 'WU713D62N9' in all_drugs? ", any(all_drugs == "WU713D62N9", na.rm=TRUE))
message("First 40 all_drugs: ", paste(head(all_drugs, 40), collapse=", "))

# 5) Read old-drug files or initialize empty
if (!is.na(olddrugfiles)) {
paths <- strsplit(olddrugfiles, ",")[[1]] %>% trimws()
old_list <- lapply(paths, function(f) {
if (file.exists(f)) fread(f, sep="\t", header=TRUE) else {
warning("Missing old‑drug file: ", f)
warning("Missing old-drug file: ", f)
NULL
}
})
Expand All @@ -101,36 +228,95 @@ if (!is.na(olddrugfiles)) {
pubchem_id=character(), canSMILES=character(),
InChIKey=character(), formula=character(), weight=numeric()
)
message("No old‑drug files provided; starting fresh")
message("No old-drug files provided; starting fresh")
}

# 6) Write placeholder
fwrite(olddrugs, newdrugfile, sep="\t", quote=FALSE)
message("Wrote placeholder to ", newdrugfile)

# 7) Augment via Python
# 7) Augment via Python (single call)
ignore_file <- "/tmp/combined_drugs_ignore_chems.txt"
use_python("/opt/venv/bin/python3", required=TRUE)
# use_python("/Users/jaco059/miniconda3/bin/python3", required=TRUE)

# source_python("build/utils/pubchem_retrieval.py")
message("\n========== DEBUG: IGNORE FILE ==========")
message("ignore_file: ", ignore_file)
message("ignore_file exists? ", file.exists(ignore_file))
if (file.exists(ignore_file)) {
ign <- tolower(trimws(readLines(ignore_file, warn=FALSE)))
message("ignore_file lines (first 50): ", paste(head(ign, 50), collapse=" | "))
message("ignore contains 'verteporfin'? ", any(ign == "verteporfin"))
message("ignore contains 'wu713d62n9'? ", any(ign == "wu713d62n9"))
}

use_python("/opt/venv/bin/python3", required=TRUE)
source_python("pubchem_retrieval.py")

# prepare prev_drug_filepaths argument for python helper (NULL if none)
prev_paths <- if (!is.na(olddrugfiles)) olddrugfiles else NULL
update_dataframe_and_write_tsv(
unique_names = all_drugs,
output_filename = newdrugfile,
ignore_chems = ignore_file,
batch_size = as.integer(50),
isname = TRUE,
prev_drug_filepaths = prev_paths,
restrict_to_raw_names = all_drugs
)


# 8) Final filter & save
tab <- fread(newdrugfile, sep="\t", header=TRUE)

# reticulate-safe list conversion
name_drugs_py <- as.list(as.character(all_drugs))

if (length(all_drugs) > 0) {
message("\n========== DEBUG: RUN SINGLE NAME-BASED PYTHON CALL ==========")
message("Calling update_dataframe_and_write_tsv with isname=TRUE for ", length(all_drugs), " names.")
message("First 30 all_drugs: ", paste(head(all_drugs, 30), collapse=", "))

update_dataframe_and_write_tsv(
unique_names = name_drugs_py,
output_filename = newdrugfile,
ignore_chems = ignore_file,
batch_size = as.integer(50),
isname = TRUE,
prev_drug_filepaths = prev_paths,
restrict_to_raw_names = name_drugs_py
)

debug_file_summary(newdrugfile, "AFTER SINGLE NAME-BASED PYTHON CALL (newdrugfile)")
} else {
message("\n========== DEBUG: SKIP PYTHON CALL (all_drugs empty) ==========")
}

# 8) Final unique & save (no CID merge)
message("\n========== DEBUG: MAIN OUTPUT BEFORE FINAL UNIQUE ==========")
debug_file_summary(newdrugfile, "MAIN OUTPUT BEFORE FINAL UNIQUE (newdrugfile)")

tab <- fread(newdrugfile, sep="\t", header=TRUE)
final_tab <- unique(tab)

message("\n========== DEBUG: FINAL CHECK ==========")
message("Final rows after unique(): ", nrow(final_tab))
if ("chem_name" %in% names(final_tab)) {
cnf <- unique(tolower(trimws(as.character(final_tab$chem_name))))
message("FINAL has verteporfin? ", any(cnf == "verteporfin", na.rm=TRUE))
message("FINAL has wu713d62n9? ", any(cnf == "wu713d62n9", na.rm=TRUE))
}

# 9) Ensure bare "verteporfin" exists as a chem_name row.
# PubChem synonyms often return only qualified variants like
# "verteporfin [usan:usp:inn:ban]" but not the bare name, which
# causes exact-join failures downstream in 03_get_experiments.
if ("chem_name" %in% names(final_tab)) {
cn_lower <- tolower(trimws(as.character(final_tab$chem_name)))
has_bare <- any(cn_lower == "verteporfin", na.rm = TRUE)

if (!has_bare) {
# Look for any row whose chem_name contains "verteporfin" or "verteporphin"
match_idx <- which(grepl("verteporfin|verteporphin", cn_lower))
if (length(match_idx) > 0) {
donor_row <- final_tab[match_idx[1], ]
donor_row[["chem_name"]] <- "verteporfin"
final_tab <- rbind(donor_row, final_tab)
message("Inserted bare 'verteporfin' alias row (from row with chem_name='",
as.character(final_tab$chem_name[match_idx[1] + 1]), # +1 because we prepended
"', improve_drug_id=", donor_row$improve_drug_id, ")")
} else {
message("WARNING: no verteporfin/verteporphin variant found in drug file; cannot insert bare alias.")
}
} else {
message("Bare 'verteporfin' already present in drug file; no insertion needed.")
}
}


fwrite(final_tab, newdrugfile, sep="\t", quote=FALSE)
message("Wrote full synonyms list to ", newdrugfile)
message("Wrote full synonyms list to ", newdrugfile)
Loading