From 1c4f112fbdba605ac323fc3a5f21cda35580a012 Mon Sep 17 00:00:00 2001 From: Jeremy Date: Mon, 21 Sep 2026 11:50:01 -0700 Subject: [PATCH 1/2] Add MPNST treated-sample and treated-omics support --- coderbuild/mpnst/00_sample_gen.R | 12 +- coderbuild/mpnst/00b_sample_gen_treated.R | 594 ++++++++++++++++++ coderbuild/mpnst/01_combined_omics.R | 32 +- coderbuild/mpnst/01b_combined_omics_treated.R | 219 +++++++ coderbuild/mpnst/build_omics.sh | 3 + coderbuild/mpnst/build_samples.sh | 2 + 6 files changed, 851 insertions(+), 11 deletions(-) create mode 100644 coderbuild/mpnst/00b_sample_gen_treated.R create mode 100644 coderbuild/mpnst/01b_combined_omics_treated.R diff --git a/coderbuild/mpnst/00_sample_gen.R b/coderbuild/mpnst/00_sample_gen.R index db1f238e..73f59c30 100644 --- a/coderbuild/mpnst/00_sample_gen.R +++ b/coderbuild/mpnst/00_sample_gen.R @@ -1,3 +1,4 @@ +# 00_sample_gen.R # This script generate a new sample table based on previous dataset's sample file (taking the max improve_sample_id) # Load required libraries library(data.table) @@ -26,9 +27,14 @@ synapser::synLogin(authToken=synapse_token) manifest<-synapser::synTableQuery("select * from syn53503360")$asDataFrame()|> as.data.frame() -#Drop contaminated sample JH-2-009 +#Drop contaminated sample JH-2-009 and others with issues manifest <- manifest %>% - filter(Sample != "JH-2-009") + filter(Sample != "JH-2-009") %>% + filter(Sample != "WU-545") %>% + filter(Sample != "WU-536") %>% + filter(Sample != "WU-505") %>% + filter(Sample != "MN-1") %>% + filter(Sample != "MN-3") ###sample file has a strict schema @@ -56,7 +62,7 @@ sampTable<-manifest|> ##third, generate a sample for the MTs if they were generated pdxmt<-subset(sampTable,!is.na(MicroTissueDrugFolder)) -pdxmt$model_type=rep('xenograft derived organoid',nrow(pdxmt)) +pdxmt$model_type=rep('3D-MEDS',nrow(pdxmt)) print(pdxmt) main<-rbind(sampTable,pdxmt)|> diff --git a/coderbuild/mpnst/00b_sample_gen_treated.R b/coderbuild/mpnst/00b_sample_gen_treated.R new file mode 100644 index 00000000..63db50b1 --- /dev/null +++ b/coderbuild/mpnst/00b_sample_gen_treated.R @@ -0,0 +1,594 @@ +# 00_sample_gen_treated.R +# This script generates treated microtissue samples and appends them to an existing samples file. + +suppressPackageStartupMessages({ + library(data.table) + library(synapser) + library(readxl) +}) + +# ----------------------- +# logging helpers +# ----------------------- +.ts <- function() format(Sys.time(), "%Y-%m-%d %H:%M:%S") +.p <- function(...) cat(.ts(), paste0(...), "\n") + +print_df_head <- function(df, n = 8) { + if (is.null(df) || nrow(df) == 0) { + .p("[PRINT] <0 rows>") + return(invisible(NULL)) + } + n <- min(n, nrow(df)) + print(utils::head(df, n)) + invisible(NULL) +} + +# ----------------------- +# constants +# ----------------------- +CONST_CANCER <- "Malignant peripheral nerve sheath tumor" +CONST_SPECIES <- "Homo sapiens (Human)" +CONST_MODEL <- "3D-MEDS" +CONST_SOURCE <- "NF Data Portal" + +# NOTE: updated to correctly preserve hyphens in IDs (MN-2, JH-2-002, WU-225) +safe_token <- function(x) { + x <- as.character(x) + x <- trimws(x) + x <- gsub("\\s+", "_", x) + x <- gsub("[^[:alnum:]_.-]+", "_", x) + x +} + +# ----------------------- +# Canonicalize individual IDs to hyphen style (MN-2, JH-2-002, WU-225) +# ----------------------- +canonicalize_individual_id <- function(x) { + x <- trimws(as.character(x)) + x <- gsub("[_\\s]+", "-", x) # underscores/spaces -> hyphen + x <- gsub("-+", "-", x) # collapse repeated hyphens + + # Normalize common patterns (only changes if it matches) + x <- sub("^MN-?([0-9A-Za-z]+)$", "MN-\\1", x) + x <- sub("^WU-?([0-9A-Za-z]+)$", "WU-\\1", x) + x <- sub("^JH-?([0-9]+)-?([0-9A-Za-z]+)$", "JH-\\1-\\2", x) + + x +} + +# ----------------------- +# Migrate existing treated other_id prefixes from underscore->hyphen +# This is because some samples were mis-labeled. +# ----------------------- +migrate_treated_other_id_prefix <- function(other_id) { + x <- as.character(other_id) + is_treated <- grepl("_treated_microtissue$", x) + + y <- x + y[is_treated] <- sub("^MN_([0-9A-Za-z]+)_", "MN-\\1_", y[is_treated]) + y[is_treated] <- sub("^WU_([0-9A-Za-z]+)_", "WU-\\1_", y[is_treated]) + y[is_treated] <- sub("^JH_([0-9]+)_([0-9A-Za-z]+)_", "JH-\\1-\\2_", y[is_treated]) + + y +} + +# ----------------------- +# Timepoint normalizer (handles "24 HR", "8HR ", etc) +# ----------------------- +normalize_timepoint <- function(tp) { + tp <- trimws(as.character(tp)) + tpu <- toupper(tp) + tpu <- gsub("\\s+", "", tpu) + + out <- tp + out[tpu %in% c("8HR", "8H")] <- "8h" + out[tpu %in% c("24HR", "24H")] <- "24h" + out[tpu %in% c("0HR", "0H")] <- "0h" + + is_h <- grepl("^[0-9]+h$", tolower(out)) + out[is_h] <- tolower(out[is_h]) + out +} + +hours_from_timepoint <- function(tp_norm) { + tp_norm <- tolower(trimws(as.character(tp_norm))) + hrs <- gsub("h$", "", tp_norm) + suppressWarnings(as.integer(hrs)) +} + +# ----------------------- +# Drop controls: DMSO + untreated (and close variants) +# ----------------------- +is_dmso <- function(drug) { + d <- tolower(trimws(as.character(drug))) + grepl("\\bdmso\\b", d) || grepl("dimethyl", d) || grepl("sulfoxide", d) +} + +is_untreated <- function(drug) { + d <- tolower(trimws(as.character(drug))) + grepl("\\buntreated\\b", d) || grepl("no[_\\s-]*treat", d) || grepl("\\bvehicle\\b", d) +} + +is_control <- function(drug) { + isTRUE(is_dmso(drug)) || isTRUE(is_untreated(drug)) +} + +# ============================================================ +# Builder function: +# - Drops DMSO + untreated rows +# - Aggregates replicates: ONE row per treated sample +# - other_id = {individual_id}_{drug}_{hours}hr_treated_microtissue +# - other_names = comma-separated replicate specimen IDs +# ============================================================ +build_from_standard_map <- function(std_map, orig_samples, batch_key, next_id, require_base = TRUE) { + std_map <- as.data.table(std_map) + + # Canonicalize IDs so other_id uses MN-2 / JH-2-002 / WU-225 formatting + std_map[, individual_id := canonicalize_individual_id(individual_id)] + + # base_exists check (optional) + std_map[, base_exists := individual_id %in% orig_samples$common_name] + + # drop controls early (DMSO + untreated) + std_map[, drop_control := vapply(drug, is_control, logical(1))] + n_ctrl <- sum(std_map$drop_control %in% TRUE) + if (n_ctrl > 0) .p("[", batch_key, "] dropping control rows (DMSO/untreated): ", n_ctrl) + std_map <- std_map[drop_control == FALSE | is.na(drop_control)] + + debug_tbl <- std_map[, .( + batch = batch_key, + specimen_id, + individual_id, + drug, + timepoint, + base_exists + )] + + keep_tbl <- if (require_base) std_map[base_exists == TRUE] else copy(std_map) + + .p("[", batch_key, "] rows after drop controls: ", nrow(std_map)) + .p("[", batch_key, "] require_base=", require_base) + .p("[", batch_key, "] rows kept pre-aggregate: ", nrow(keep_tbl)) + .p("[", batch_key, "] preview:") + print_df_head(keep_tbl, 8) + + if (nrow(keep_tbl) == 0) { + return(list( + new_samples_tbl = data.table(), + map_tbl = data.table(), + debug_tbl = debug_tbl, + next_id = next_id + )) + } + + # ---- First aggregation: collapse exact (individual_id, drug, timepoint) rows + agg1 <- keep_tbl[ + , + .(specimen_ids = sort(unique(trimws(as.character(specimen_id))))), + by = .(individual_id, drug, timepoint) + ] + + # Remove NA/blank specimen IDs inside groups + agg1[, specimen_ids := lapply(specimen_ids, function(v) v[!is.na(v) & nzchar(v)])] + agg1 <- agg1[lengths(specimen_ids) > 0] + + agg1[, other_names := vapply(specimen_ids, function(v) paste(v, collapse = ", "), character(1))] + agg1[, hours := hours_from_timepoint(timepoint)] + agg1[is.na(hours), hours := -1L] + + # ---- Compute other_id_val per group (this is what ends up in samples.csv) + agg1[, other_id_val := paste0( + safe_token(individual_id), "_", + safe_token(drug), "_", + as.integer(hours), "hr_treated_microtissue" + )] + + # ---- SECOND aggregation: merge replicate IDs when other_id_val collides + agg <- agg1[ + , + .( + individual_id = individual_id[1], + drug = drug[1], + timepoint = timepoint[1], + hours = hours[1], + other_names = paste( + sort(unique(trimws(unlist(strsplit(paste(other_names, collapse = ", "), ",\\s*"))))), + collapse = ", " + ) + ), + by = .(other_id_val) + ] + + .p("[", batch_key, "] aggregate groups after merge-by-other_id: ", nrow(agg)) + .p("[", batch_key, "] aggregate head:") + print_df_head(agg[, .(other_id_val, individual_id, drug, timepoint, other_names, hours)], 8) + + # fast lookup of existing treated keys in orig_samples + orig_treated_keys <- paste0(as.character(orig_samples$other_id), "||", as.character(orig_samples$model_type)) + + map_tbl <- data.table( + batch = character(), + individual_id = character(), + drug = character(), + timepoint = character(), + hours = integer(), + other_id = character(), + other_names = character(), + treated_improve_sample_id = integer() + ) + + new_samples_tbl <- data.table( + other_id = character(), + common_name = character(), + other_id_source = character(), + other_names = character(), + cancer_type = character(), + species = character(), + model_type = character(), + improve_sample_id = integer() + ) + + for (i in seq_len(nrow(agg))) { + indiv <- as.character(agg$individual_id[i]) + drug <- as.character(agg$drug[i]) + tp <- as.character(agg$timepoint[i]) + hrs <- as.integer(agg$hours[i]) + repnames <- as.character(agg$other_names[i]) + + other_id_val <- as.character(agg$other_id_val[i]) + key <- paste0(other_id_val, "||", CONST_MODEL) + + # Skip if already exists in samples.csv + if (key %in% orig_treated_keys) { + .p("[", batch_key, "][SKIP] already exists other_id=", other_id_val) + next + } + + treated_id <- next_id + next_id <- next_id + 1L + + new_samples_tbl <- rbind( + new_samples_tbl, + data.table( + other_id = other_id_val, + common_name = indiv, + other_id_source = CONST_SOURCE, + other_names = repnames, + cancer_type = CONST_CANCER, + species = CONST_SPECIES, + model_type = CONST_MODEL, + improve_sample_id = treated_id + ), + fill = TRUE + ) + + map_tbl <- rbind( + map_tbl, + data.table( + batch = batch_key, + individual_id = indiv, + drug = drug, + timepoint = tp, + hours = hrs, + other_id = other_id_val, + other_names = repnames, + treated_improve_sample_id = treated_id + ), + fill = TRUE + ) + + .p("[", batch_key, "][ADD] ", other_id_val, " | reps={", repnames, "} -> ", treated_id) + } + + list( + new_samples_tbl = new_samples_tbl, + map_tbl = map_tbl, + debug_tbl = debug_tbl, + next_id = next_id + ) +} + +write_batch_outputs <- function(batch_key, debug_tbl, map_tbl) { + out_debug_path <- file.path("/tmp", sprintf("mpnst_samples_treated_%s_debug.csv", batch_key)) + out_map_path <- file.path("/tmp", sprintf("mpnst_samples_treated_%s_map.csv", batch_key)) + + .p("[WRITE] ", batch_key, " debug -> ", out_debug_path) + fwrite(debug_tbl, out_debug_path) + .p("[WRITE] OK") + + .p("[WRITE] ", batch_key, " map -> ", out_map_path) + fwrite(map_tbl, out_map_path) + .p("[WRITE] OK") +} + +# ============================================================ +# Batch 1 +# ============================================================ +BATCH1_KEY <- "batch1" +BATCH1_SAMPLEMAP_SYNID <- "syn64608368" + +fetch_batch1_map <- function() { + .p("[BATCH1] synGet(", BATCH1_SAMPLEMAP_SYNID, ") ...") + pth <- synGet(BATCH1_SAMPLEMAP_SYNID)$path + dt <- as.data.table(fread(pth)) + .p("[BATCH1] map dim: ", nrow(dt), " x ", ncol(dt)) + print_df_head(dt, 6) + dt +} + +normalize_batch1_map <- function(raw_map) { + dt <- copy(as.data.table(raw_map)) + dt[, specimenID := trimws(as.character(specimenID))] + dt[, individualID := trimws(as.character(individualID))] + dt[, Drug := trimws(as.character(Drug))] + dt[, Timepoint := trimws(as.character(Timepoint))] + + dt <- dt[!is.na(specimenID) & nzchar(specimenID) & !is.na(individualID) & nzchar(individualID)] + + dt[, .( + specimen_id = specimenID, + individual_id = individualID, + drug = Drug, + timepoint = normalize_timepoint(Timepoint) + )] +} + +run_batch1 <- function(orig_samples, next_id) { + build_from_standard_map(normalize_batch1_map(fetch_batch1_map()), orig_samples, BATCH1_KEY, next_id, require_base = TRUE) +} + +# ============================================================ +# Batch 2 +# ============================================================ +BATCH2_KEY <- "batch2" +BATCH2_SAMPLEMAP_SYNID <- "syn66302373" # Update to syn64608372 once the corrected version is uploaded. + + +fetch_batch2_map <- function() { + .p("[BATCH2] synGet(", BATCH2_SAMPLEMAP_SYNID, ") ...") + pth <- synGet(BATCH2_SAMPLEMAP_SYNID)$path + dt <- as.data.table(fread(pth)) + .p("[BATCH2] map dim: ", nrow(dt), " x ", ncol(dt)) + print_df_head(dt, 6) + dt +} + +normalize_batch2_map <- function(raw_map) { + dt <- copy(as.data.table(raw_map)) + dt[, specimenID := trimws(as.character(specimenID))] + dt[, individualID := trimws(as.character(individualID))] + dt[, Drug := trimws(as.character(Drug))] + dt[, Timepoint := trimws(as.character(Timepoint))] + + #fixed in metadata, so this is not needed anymore. + # MN2 -> MN-2 + # dt[individualID == "MN2", individualID := "MN-2"] + + dt <- dt[!is.na(specimenID) & nzchar(specimenID) & !is.na(individualID) & nzchar(individualID)] + + dt[, .( + specimen_id = specimenID, + individual_id = individualID, + drug = Drug, + timepoint = normalize_timepoint(Timepoint) + )] +} + +run_batch2 <- function(orig_samples, next_id) { + build_from_standard_map(normalize_batch2_map(fetch_batch2_map()), orig_samples, BATCH2_KEY, next_id, require_base = TRUE) +} + +# ============================================================ +# Batch 3 +# ============================================================ +BATCH3_KEY <- "batch3" +BATCH3_SAMPLEMAP_SYNID <- "syn66050299" + +fetch_batch3_map <- function() { + .p("[BATCH3] synGet(", BATCH3_SAMPLEMAP_SYNID, ") ...") + pth <- synGet(BATCH3_SAMPLEMAP_SYNID)$path + dt <- as.data.table(fread(pth)) + .p("[BATCH3] map dim: ", nrow(dt), " x ", ncol(dt)) + print_df_head(dt, 6) + dt +} + +normalize_batch3_map <- function(raw_map) { + dt <- copy(as.data.table(raw_map)) + dt[, specimenID := trimws(as.character(specimenID))] + dt[, individualID := trimws(as.character(individualID))] + dt[, Drug := trimws(as.character(Drug))] + dt[, Timepoint := trimws(as.character(Timepoint))] + + dt <- dt[!is.na(specimenID) & nzchar(specimenID) & !is.na(individualID) & nzchar(individualID)] + + dt[, .( + specimen_id = specimenID, + individual_id = individualID, + drug = Drug, + timepoint = normalize_timepoint(Timepoint) + )] +} + +run_batch3 <- function(orig_samples, next_id) { + build_from_standard_map(normalize_batch3_map(fetch_batch3_map()), orig_samples, BATCH3_KEY, next_id, require_base = TRUE) +} + +# ============================================================ +# Batch 4 (CSV) +# ============================================================ +BATCH4_KEY <- "batch4" +BATCH4_SAMPLEMAP_SYNID <- "syn72518652" + +fetch_batch4_map <- function() { + .p("[BATCH4] synGet(", BATCH4_SAMPLEMAP_SYNID, ") ...") + pth <- synGet(BATCH4_SAMPLEMAP_SYNID)$path + dt <- as.data.table(fread(pth)) + .p("[BATCH4] map dim: ", nrow(dt), " x ", ncol(dt)) + print_df_head(dt, 6) + dt +} + +normalize_batch4_map <- function(raw_map) { + dt <- copy(as.data.table(raw_map)) + + # Expected columns in 01252026_Batch4_samplemap.csv: + # specimenID, individualID, Drug, Timepoint + if (!all(c("specimenID", "individualID", "Drug", "Timepoint") %in% names(dt))) { + .p("[BATCH4][WARN] Missing required columns. Found: ", paste(names(dt), collapse = " | ")) + return(data.table( + specimen_id = character(), + individual_id = character(), + drug = character(), + timepoint = character() + )) + } + + dt[, specimenID := trimws(as.character(specimenID))] + dt[, individualID := trimws(as.character(individualID))] + dt[, Drug := trimws(as.character(Drug))] + dt[, Timepoint := trimws(as.character(Timepoint))] + + dt <- dt[!is.na(specimenID) & nzchar(specimenID) & !is.na(individualID) & nzchar(individualID)] + + out <- dt[, .( + specimen_id = specimenID, + individual_id = individualID, + drug = Drug, + timepoint = normalize_timepoint(Timepoint) + )] + + .p("[BATCH4] normalized dim: ", nrow(out), " x ", ncol(out)) + .p("[BATCH4] normalized head:") + print_df_head(out, 8) + + out +} + +run_batch4 <- function(orig_samples, next_id) { + build_from_standard_map( + normalize_batch4_map(fetch_batch4_map()), + orig_samples, + BATCH4_KEY, + next_id, + require_base = FALSE + ) +} + + +# ============================================================ +# Register all batches (order) +# ============================================================ +ALL_BATCH_RUNNERS <- list( + list(key=BATCH1_KEY, fn=run_batch1), + list(key=BATCH2_KEY, fn=run_batch2), + list(key=BATCH3_KEY, fn=run_batch3), + list(key=BATCH4_KEY, fn=run_batch4) +) + +# ============================================================ +# MAIN +# ============================================================ +args <- commandArgs(trailingOnly = TRUE) + +.p("============================================================") +.p("=== 00b_sample_gen_treated.R (RUN ALL BATCHES) ===") +.p("Timestamp: ", .ts()) +.p("Args: ", paste(args, collapse = " | ")) +.p("============================================================") + +if (length(args) < 1) { + stop( + "Usage: Rscript 00b_sample_gen_treated.R \n", + "Example: Rscript 00b_sample_gen_treated.R /tmp/mpnst_samples.csv", + call. = FALSE + ) +} + +prev_samples_path <- args[1] +.p("[ARGS] prev_samples_path=", prev_samples_path) +if (!file.exists(prev_samples_path)) stop("File not found: ", prev_samples_path, call. = FALSE) + +.p("[READ] fread(samples) ...") +orig_samples <- fread(prev_samples_path) +.p("[READ] orig_samples dim: ", nrow(orig_samples), " x ", ncol(orig_samples)) +print_df_head(orig_samples, 6) + +# ------------------------------------------------------------ +# Migrate already-existing treated other_ids to hyphen style +# so future runs match your desired format and don't skip forever +# ------------------------------------------------------------ +orig_samples[, other_id_old := other_id] +orig_samples[, other_id := migrate_treated_other_id_prefix(other_id)] + +n_changed <- sum(orig_samples$other_id != orig_samples$other_id_old, na.rm = TRUE) +if (n_changed > 0) .p("[MIGRATE] updated treated other_id prefix for ", n_changed, " row(s)") +orig_samples[, other_id_old := NULL] + +# If migration created duplicates (same other_id + model_type), merge other_names and drop dup rows +orig_samples[, dedup_key := paste(other_id, model_type, sep = "||")] +dup_keys <- orig_samples[duplicated(dedup_key), unique(dedup_key)] + +if (length(dup_keys) > 0) { + .p("[MIGRATE] duplicates after migration: ", length(dup_keys), " key(s). Merging other_names + dropping dup rows.") + + merged <- orig_samples[dedup_key %in% dup_keys, .( + merged_other_names = paste( + sort(unique(trimws(unlist(strsplit(paste(na.omit(other_names), collapse = ", "), ",\\s*"))))), + collapse = ", " + ) + ), by = dedup_key] + + orig_samples[merged, on = "dedup_key", other_names := i.merged_other_names] + orig_samples <- orig_samples[!duplicated(dedup_key)] +} + +orig_samples[, dedup_key := NULL] + +token <- Sys.getenv("SYNAPSE_AUTH_TOKEN") +.p("[ENV] SYNAPSE_AUTH_TOKEN present? ", nzchar(token)) +if (!nzchar(token)) stop("Please set SYNAPSE_AUTH_TOKEN in your environment.", call. = FALSE) + +.p("[SYN] synLogin() ...") +synLogin(authToken = token) +.p("[SYN] synLogin() OK") + +max_id <- suppressWarnings(max(as.integer(orig_samples$improve_sample_id), na.rm = TRUE)) +if (!is.finite(max_id)) max_id <- 0L +next_id <- max_id + 1L +.p("[IDS] max improve_sample_id=", max_id, " next_id=", next_id) + +.p("============================================================") +.p("[RUN] Running ", length(ALL_BATCH_RUNNERS), " batch runner(s) in order...") +.p("============================================================") + +for (item in ALL_BATCH_RUNNERS) { + batch_key <- item$key + runner <- item$fn + + res <- runner(orig_samples = orig_samples, next_id = next_id) + + if (!is.null(res$new_samples_tbl) && nrow(res$new_samples_tbl) > 0) { + orig_samples <- rbindlist(list(orig_samples, res$new_samples_tbl), use.names = TRUE, fill = TRUE) + + orig_samples[, dedup_key := paste(other_id, model_type, sep = "||")] + before <- nrow(orig_samples) + orig_samples <- orig_samples[!duplicated(dedup_key)] + after <- nrow(orig_samples) + if (after < before) .p("[", batch_key, "][WARN] post-append dedup removed ", before - after, " row(s)") + orig_samples[, dedup_key := NULL] + } + + next_id <- res$next_id + + .p("[RUN] Finished ", batch_key, " | added=", ifelse(is.null(res$new_samples_tbl), 0, nrow(res$new_samples_tbl))) + write_batch_outputs(batch_key, res$debug_tbl, res$map_tbl) +} + +.p("[WRITE] Writing updated samples -> ", prev_samples_path) +fwrite(orig_samples, prev_samples_path) +.p("[WRITE] OK") + +.p("============================================================") +.p("[DONE] Completed all batches. Updated samples at: ", prev_samples_path) +.p("============================================================") diff --git a/coderbuild/mpnst/01_combined_omics.R b/coderbuild/mpnst/01_combined_omics.R index dcbdfbae..7e7dd41a 100644 --- a/coderbuild/mpnst/01_combined_omics.R +++ b/coderbuild/mpnst/01_combined_omics.R @@ -1,3 +1,4 @@ +# 01_combined_omics.R #!/usr/bin/env Rscript # Combined MPNST & MPNST-PDX Data Extraction Script @@ -31,7 +32,7 @@ genes_df <- fread(genes) # Subset by model type pdx_samps <- filter(samples_df, model_type == "patient derived xenograft") tumor_samps<- filter(samples_df, model_type == "tumor") -mt_samps <- filter(samples_df, model_type == "xenograft derived organoid") # These end up being the same as pdx_samps in the manifest. +mt_samps <- filter(samples_df, model_type == "3D-MEDS") # These end up being the same as pdx_samps in the manifest. # Retrieve manifest table from Synapse manifest <- synTableQuery("select * from syn53503360")$asDataFrame() %>% @@ -58,7 +59,7 @@ tumor_data <- manifest %>% mutate(Proteomics = "") %>% filter(!is.na(improve_sample_id)) -mt_data <- manifest %>% #Note, this is the same as pdx_data but I think we default to "xenograft derived organoid" if present (based on original files) +mt_data <- manifest %>% #Note, this is the same as pdx_data but I think we default to "3D-MEDS" if present (based on original files) select(common_name, starts_with("PDX")) %>% left_join(mt_samps, by = "common_name") %>% select(improve_sample_id, common_name, model_type, @@ -78,7 +79,7 @@ study_label <- function(type) { case_when( type == "patient derived xenograft" ~ "MPNST PDX", type == "tumor" ~ "MPNST Tumor", - type == "xenograft derived organoid" ~ "MPNST PDX MT", + type == "3D-MEDS" ~ "MPNST PDX MT", TRUE ~ "MPNST" ) } @@ -147,13 +148,28 @@ transcriptomics_list <- lapply( if (is.null(meta)) return(NULL) df <- tryCatch({ - fread(synGet(id)$path) %>% + raw <- fread(synGet(id)$path) %>% separate(Name, into = c("other_id","vers"), sep = "\\.") %>% - select(-vers) %>% + select(-vers) + + is_enst_input <- all(grepl("^ENST", raw$other_id)) + + mapped <- raw %>% left_join(genes_df) %>% - select(entrez_id, transcriptomics = TPM) %>% - filter(!is.na(entrez_id), transcriptomics != 0) %>% - distinct() + select(entrez_id, other_id, transcriptomics = TPM) %>% + filter(!is.na(entrez_id), transcriptomics != 0) + + if (is_enst_input) { + mapped <- mapped %>% + group_by(entrez_id) %>% + summarise(transcriptomics = sum(transcriptomics), .groups = "drop") + } else { + mapped <- mapped %>% + select(entrez_id, transcriptomics) %>% + distinct() + } + + mapped }, error = function(e) NULL) i_safe_extract( diff --git a/coderbuild/mpnst/01b_combined_omics_treated.R b/coderbuild/mpnst/01b_combined_omics_treated.R new file mode 100644 index 00000000..8125fc61 --- /dev/null +++ b/coderbuild/mpnst/01b_combined_omics_treated.R @@ -0,0 +1,219 @@ +# 01b_combined_omics_treated.R +#!/usr/bin/env Rscript + +suppressPackageStartupMessages({ + library(data.table) + library(synapser) + library(dplyr) + library(tidyr) +}) + +# ----------------------- +# logging helpers +# ----------------------- +.ts <- function() format(Sys.time(), "%Y-%m-%d %H:%M:%S") +.p <- function(...) cat(.ts(), paste0(...), "\n") + +# ----------------------- +# batch inputs (Synapse) +# ----------------------- +B1_SYN <- "syn64608366" # salmon.merged.gene_tpm_corrected_id.csv +B2_SYN <- "syn64608370" # salmon.merged.gene_tpm.tsv_corrected_ID.csv +B3_SYN <- "syn65887795" # salmon.merged.gene_tpm.tsv +B4_SYN <- "syn68898091" # salmon.merged.gene_tpm.tsv + +# ----------------------- +# helpers +# ----------------------- +study_label <- function(type) { + dplyr::case_when( + type == "patient derived xenograft" ~ "MPNST PDX", + type == "tumor" ~ "MPNST Tumor", + type == "3D-MEDS" ~ "MPNST PDX MT", + TRUE ~ "MPNST" + ) +} + +is_treated_microtissue <- function(other_id, other_names) { + oid <- as.character(other_id) + onm <- as.character(other_names) + grepl("treated_microtissue$", oid) && !is.na(onm) && nzchar(trimws(onm)) +} + +split_rep_ids <- function(x) { + if (is.na(x) || !nzchar(trimws(x))) return(character(0)) + v <- unlist(strsplit(x, ",")) + v <- trimws(v) + v <- v[nzchar(v)] + unique(v) +} + +guess_batch_from_rep <- function(rep_id) { + r <- as.character(rep_id) + if (grepl("^b1_", r)) return("b1") + if (grepl("^b2_", r)) return("b2") + if (grepl("^B3-", r) || grepl("^B3_", r)) return("b3") + if (grepl("^B4-", r) || grepl("^B4\\.", r)) return("b4") + return(NA_character_) +} + +rep_to_colname <- function(rep_id, batch) { + r <- as.character(rep_id) + if (batch %in% c("b1", "b2")) return(r) + if (batch == "b3") return(gsub("-", "_", r, fixed = FALSE)) # B3-10 -> B3_10 + if (batch == "b4") return(gsub("B4-", "B4.", r, fixed = TRUE)) # B4-10 -> B4.10 + return(r) +} + +read_tpm_matrix <- function(syn_id) { + .p("[SYN] synGet(", syn_id, ") ...") + pth <- synGet(syn_id)$path + .p("[READ] ", syn_id, " -> ", pth) + dt <- fread(pth) + stopifnot(all(c("gene_id", "gene_name") %in% names(dt))) + dt +} + +# returns a data.table with columns: +# gene_id, transcriptomics, improve_sample_id, source, study +extract_treated_from_batch <- function(mat, treated_tbl, batch_key) { + if (nrow(treated_tbl) == 0) return(data.table()) + + tmp <- copy(treated_tbl) + tmp[, rep_ids := lapply(other_names, split_rep_ids)] + tmp[, rep_batch := vapply(rep_ids, function(v) { + if (length(v) == 0) return(NA_character_) + guess_batch_from_rep(v[1]) + }, character(1))] + + tmp <- tmp[rep_batch == batch_key] + if (nrow(tmp) == 0) return(data.table()) + + tmp[, rep_cols := lapply(rep_ids, function(v) rep_to_colname(v, batch_key))] + + .p("[", batch_key, "] treated samples in this batch: ", nrow(tmp)) + + out_list <- vector("list", nrow(tmp)) + + for (i in seq_len(nrow(tmp))) { + samp_id <- tmp$improve_sample_id[i] + samp_model <- tmp$model_type[i] + cols <- tmp$rep_cols[[i]] + cols <- cols[cols %in% names(mat)] + + if (length(cols) == 0) { + .p("[", batch_key, "][SKIP] improve_sample_id=", samp_id, + " -> no replicate columns found in matrix") + next + } + + expr <- rowMeans(as.matrix(mat[, ..cols]), na.rm = TRUE) + + out_list[[i]] <- data.table( + gene_id = mat$gene_id, + transcriptomics = expr, + improve_sample_id = samp_id, + source = "NF Data Portal", + study = study_label(samp_model) + ) + } + + rbindlist(out_list, use.names = TRUE, fill = TRUE) +} + +# ----------------------- +# MAIN +# ----------------------- +args <- commandArgs(trailingOnly = TRUE) +if (length(args) < 3) { + stop("Usage: Rscript 01b_combined_omics_treated.R ", call. = FALSE) +} +PAT <- args[1] +samples <- args[2] +genes <- args[3] + +.p("============================================================") +.p("=== 01b_combined_omics_treated.R ===") +.p("Timestamp: ", .ts()) +.p("Args: ", paste(args, collapse = " | ")) +.p("============================================================") + +.p("[SYN] synLogin() ...") +synLogin(authToken = PAT) +.p("[SYN] synLogin() OK") + +.p("[READ] samples: ", samples) +samples_df <- fread(samples) + +treated <- samples_df[ + model_type == "3D-MEDS" & + vapply(seq_len(nrow(samples_df)), function(i) { + is_treated_microtissue(samples_df$other_id[i], samples_df$other_names[i]) + }, logical(1)), + .(other_names, improve_sample_id, model_type) +] + +.p("[SAMPLES] treated_microtissue rows: ", nrow(treated)) +if (nrow(treated) == 0) stop("No treated microtissue samples found in samples.csv", call. = FALSE) + +.p("[READ] genes: ", genes) +genes_df <- fread(genes) +stopifnot(all(c("other_id", "entrez_id") %in% names(genes_df))) + +# matrices +m1 <- read_tpm_matrix(B1_SYN) +m2 <- read_tpm_matrix(B2_SYN) +m3 <- read_tpm_matrix(B3_SYN) +m4 <- read_tpm_matrix(B4_SYN) + +# extract +t1 <- extract_treated_from_batch(m1, treated, "b1") +t2 <- extract_treated_from_batch(m2, treated, "b2") +t3 <- extract_treated_from_batch(m3, treated, "b3") +t4 <- extract_treated_from_batch(m4, treated, "b4") + +treated_long <- rbindlist(list(t1, t2, t3, t4), use.names = TRUE, fill = TRUE) +.p("[OUT] treated_long rows: ", nrow(treated_long)) +if (nrow(treated_long) == 0) stop("No treated data extracted (check replicate IDs vs column names).", call. = FALSE) + +# ENSEMBL gene_id has version -> strip to match genes_df$other_id +treated_long[, other_id := tstrsplit(gsub("\\s+", "", gene_id), "\\.", keep = 1)] + +treated_final <- treated_long %>% + left_join(genes_df, by = "other_id") %>% + dplyr::select(entrez_id, transcriptomics, improve_sample_id, source, study) %>% + filter(!is.na(entrez_id), transcriptomics != 0) %>% + distinct() + +# enforce strict schema & order +treated_final <- as.data.table(treated_final)[, .(entrez_id, transcriptomics, improve_sample_id, source, study)] + +treated_out_path <- file.path("/tmp", "mpnst_transcriptomics_treated_microtissue.csv") +.p("[WRITE] treated-only transcriptomics -> ", treated_out_path) +fwrite(treated_final, treated_out_path) +.p("[WRITE] OK") + +# optional append into main transcriptomics (same strict schema) +main_path <- file.path("/tmp", "mpnst_transcriptomics.csv") +if (file.exists(main_path)) { + .p("[APPEND] Found existing main transcriptomics: ", main_path) + main_df <- fread(main_path) + main_df <- main_df[, .(entrez_id, transcriptomics, improve_sample_id, source, study)] # enforce schema + + combined <- rbindlist(list(main_df, treated_final), use.names = TRUE, fill = TRUE) + + # de-dup conservatively on (entrez_id, improve_sample_id, source, study, transcriptomics) + combined[, dedup_key := paste(entrez_id, improve_sample_id, source, study, transcriptomics, sep = "||")] + combined <- combined[!duplicated(dedup_key)] + combined[, dedup_key := NULL] + + .p("[WRITE] overwriting main transcriptomics with appended treated rows -> ", main_path) + fwrite(combined, main_path) + .p("[WRITE] OK") +} else { + .p("[APPEND] No existing /tmp/mpnst_transcriptomics.csv found. Not appending.") +} + +.p("============================================================") +.p("[DONE] Wrote treated transcriptomics with strict schema.") +.p("============================================================") diff --git a/coderbuild/mpnst/build_omics.sh b/coderbuild/mpnst/build_omics.sh index d6d2cec7..c6b0853c 100644 --- a/coderbuild/mpnst/build_omics.sh +++ b/coderbuild/mpnst/build_omics.sh @@ -1,3 +1,4 @@ +# build_omics.sh #!/bin/bash set -euo pipefail @@ -5,3 +6,5 @@ trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit echo "Running 01_combined_omics.R with $SYNAPSE_AUTH_TOKEN, $2, and $1." Rscript 01_combined_omics.R $SYNAPSE_AUTH_TOKEN $2 $1 +echo "Running 01b_combined_omics_treated.R with PAT, samples=$2, genes=$1" +Rscript 01b_combined_omics_treated.R "$SYNAPSE_AUTH_TOKEN" "$2" "$1" \ No newline at end of file diff --git a/coderbuild/mpnst/build_samples.sh b/coderbuild/mpnst/build_samples.sh index c9c079fa..832ed895 100644 --- a/coderbuild/mpnst/build_samples.sh +++ b/coderbuild/mpnst/build_samples.sh @@ -1,3 +1,4 @@ +# build_samples.sh #!/bin/bash set -euo pipefail @@ -5,3 +6,4 @@ trap 'echo "Error on or near line $LINENO while executing: $BASH_COMMAND"; exit echo "Running 00_sample_gen.R with $1." Rscript 00_sample_gen.R $1 +Rscript 00b_sample_gen_treated.R /tmp/mpnst_samples.csv From 55ec7e167c8b0750d1c1ed4ecbae1bc7ea01fcab Mon Sep 17 00:00:00 2001 From: Jeremy Date: Mon, 21 Sep 2026 11:50:01 -0700 Subject: [PATCH 2/2] Add the cNF organoid dataset: drugs and experiments --- coderbuild/cnf/03-drugs-cnf.py | 531 +++++++++++++++++++++++++++ coderbuild/cnf/04-experiments-cnf.py | 478 ++++++++++++++++++++++++ coderbuild/cnf/build_drugs.sh | 20 + coderbuild/cnf/build_exp.sh | 12 + 4 files changed, 1041 insertions(+) create mode 100644 coderbuild/cnf/03-drugs-cnf.py create mode 100644 coderbuild/cnf/04-experiments-cnf.py create mode 100644 coderbuild/cnf/build_drugs.sh create mode 100644 coderbuild/cnf/build_exp.sh diff --git a/coderbuild/cnf/03-drugs-cnf.py b/coderbuild/cnf/03-drugs-cnf.py new file mode 100644 index 00000000..fe307ddb --- /dev/null +++ b/coderbuild/cnf/03-drugs-cnf.py @@ -0,0 +1,531 @@ +#!/usr/bin/env python3 +""" +03-drugs-cnf.py: generate cnf_drugs.tsv and cnf_drug_descriptors.tsv. + +Pulls unique drug names from the cNF drug screen on Synapse, then calls the +standard coderdata utilities: + + coderbuild/utils/pubchem_retrieval.py + -> cnf_drugs.tsv + + coderbuild/utils/build_drug_descriptor_table.py + -> cnf_drug_descriptors.tsv + +Important implementation note: + pubchem_retrieval.py does not expose a command-line interface in the + provided version. It defines update_dataframe_and_write_tsv(), but running + `python pubchem_retrieval.py --input ... --output ...` exits 0 without + doing any work. Therefore this script imports pubchem_retrieval.py directly + and calls update_dataframe_and_write_tsv(). + +The descriptor builder is still run as a subprocess. It may call R/biomaRt +internally, which can fail when Ensembl is unavailable. To keep the drugs file +produced even when descriptors fail, this script: + + 1. Retries external/remote work with exponential backoff. + 2. Sets a generous timeout. + 3. Supports --skip-descriptors. + 4. Supports --only-descriptors for retrying descriptor generation later. + +Usage: + python 03-drugs-cnf.py --prev_drugs + [--out_drugs /tmp/cnf_drugs.tsv] + [--out_desc /tmp/cnf_drug_descriptors.tsv] + [--skip-descriptors] + [--only-descriptors] + [--retries 3] + [--timeout 1800] +""" + +import argparse +import importlib.util +import logging +import os +import signal +import subprocess +import sys +import tempfile +import time +import pandas as pd +import synapseclient + + +# ---------------------------------------------------------------------------- +# Constants +# ---------------------------------------------------------------------------- + +DRUG_SCREEN_TABLE = "syn51301431" + +PUBCHEM_SCRIPT = os.environ.get( + "CODERDATA_PUBCHEM_SCRIPT", + "/coderbuild/utils/pubchem_retrieval.py", +) + +DESCRIPTOR_SCRIPT = os.environ.get( + "CODERDATA_DESCRIPTOR_SCRIPT", + "/coderbuild/utils/build_drug_desc.py", +) + +DEFAULT_OUT_DRUGS = "/tmp/cnf_drugs.tsv" +DEFAULT_OUT_DESC = "/tmp/cnf_drug_descriptors.tsv" + +DEFAULT_RETRIES = 3 +DEFAULT_TIMEOUT = 1800 # seconds per attempt +RETRY_BACKOFF_S = 30 # 30s, 60s, 120s, ... + + +# ---------------------------------------------------------------------------- +# Logging +# ---------------------------------------------------------------------------- + +def configure_logging() -> None: + logging.basicConfig( + format="[%(asctime)s] [%(levelname)s] %(message)s", + datefmt="%Y-%m-%d %H:%M:%S", + level=logging.INFO, + ) + + +# ---------------------------------------------------------------------------- +# Generic subprocess wrapper with retry +# ---------------------------------------------------------------------------- + +def run_with_retries( + cmd: list[str], + retries: int, + timeout: int, + label: str, +) -> bool: + """ + Run a subprocess, retrying on failure with exponential backoff. + + Returns True on success, False if all retries are exhausted. + """ + for attempt in range(1, retries + 1): + logging.info("[%s] attempt %d/%d", label, attempt, retries) + + try: + subprocess.run(cmd, check=True, timeout=timeout) + logging.info("[%s] succeeded on attempt %d", label, attempt) + return True + + except subprocess.TimeoutExpired: + logging.warning( + "[%s] timed out after %ds on attempt %d", + label, + timeout, + attempt, + ) + + except subprocess.CalledProcessError as exc: + logging.warning( + "[%s] failed with exit %d on attempt %d", + label, + exc.returncode, + attempt, + ) + + except Exception as exc: + logging.warning( + "[%s] unexpected error on attempt %d: %s", + label, + attempt, + exc, + ) + + if attempt < retries: + backoff = RETRY_BACKOFF_S * (2 ** (attempt - 1)) + logging.info("[%s] sleeping %ds before retry", label, backoff) + time.sleep(backoff) + + logging.error("[%s] exhausted %d retries", label, retries) + return False + + +# ---------------------------------------------------------------------------- +# PubChem utility loading/calling +# ---------------------------------------------------------------------------- + +def load_pubchem_module(path: str): + """ + Load pubchem_retrieval.py as a Python module. + + The provided pubchem_retrieval.py has no command-line entry point, so it + must be imported and called directly. + """ + if not os.path.exists(path): + raise FileNotFoundError(f"PubChem retrieval script not found: {path}") + + spec = importlib.util.spec_from_file_location( + "coderdata_pubchem_retrieval", + path, + ) + + if spec is None or spec.loader is None: + raise ImportError(f"Could not create import spec for {path}") + + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + + if not hasattr(module, "update_dataframe_and_write_tsv"): + raise AttributeError( + f"{path} does not define update_dataframe_and_write_tsv()" + ) + + return module + + +def call_pubchem_retrieval( + names: list[str], + prev_drug_files: list[str], + out_drugs: str, + retries: int, + timeout: int, +) -> bool: + """ + Generate the drugs TSV by calling pubchem_retrieval.update_dataframe_and_write_tsv(). + + This intentionally does not use subprocess because the supplied + pubchem_retrieval.py has no argparse/main block and exits 0 without writing + anything when run as a script. + """ + prev_arg = ",".join(prev_drug_files) if prev_drug_files else None + + for attempt in range(1, retries + 1): + logging.info("[pubchem_retrieval] attempt %d/%d", attempt, retries) + + ignore_chems = None + + try: + pubchem = load_pubchem_module(PUBCHEM_SCRIPT) + + with tempfile.NamedTemporaryFile( + "w", + suffix="_cnf_ignore_chems.txt", + delete=False, + ) as tmp_ignore: + ignore_chems = tmp_ignore.name + + df = pubchem.update_dataframe_and_write_tsv( + unique_names=names, + output_filename=out_drugs, + ignore_chems=ignore_chems, + batch_size=1, + isname=True, + time_limit=timeout, + prev_drug_filepaths=prev_arg, + restrict_to_raw_names=names, + ) + + # pubchem_retrieval.py sets signal.alarm(time_limit), but does not + # clear it. Clear it here so the descriptor step is not interrupted. + try: + signal.alarm(0) + except Exception: + pass + + if getattr(pubchem, "should_continue", True) is False: + raise TimeoutError( + f"pubchem_retrieval reached its time limit of {timeout}s" + ) + + if not os.path.exists(out_drugs): + raise FileNotFoundError( + f"pubchem_retrieval returned but did not create {out_drugs}" + ) + + if os.path.getsize(out_drugs) == 0: + raise RuntimeError( + f"pubchem_retrieval created {out_drugs}, but it is empty" + ) + + n_rows = len(df) if df is not None else "unknown" + + logging.info( + "[pubchem_retrieval] succeeded on attempt %d; wrote %s rows to %s", + attempt, + n_rows, + out_drugs, + ) + + return True + + except Exception as exc: + try: + signal.alarm(0) + except Exception: + pass + + logging.warning( + "[pubchem_retrieval] failed on attempt %d/%d: %s", + attempt, + retries, + exc, + ) + + if attempt < retries: + backoff = RETRY_BACKOFF_S * (2 ** (attempt - 1)) + logging.info( + "[pubchem_retrieval] sleeping %ds before retry", + backoff, + ) + time.sleep(backoff) + + finally: + if ignore_chems and os.path.exists(ignore_chems): + try: + os.unlink(ignore_chems) + except OSError: + pass + + logging.error("[pubchem_retrieval] exhausted %d retries", retries) + return False + + +# ---------------------------------------------------------------------------- +# Drug name collection +# ---------------------------------------------------------------------------- + +def collect_drug_names(syn) -> list[str]: + """ + Pull every cNF drug-screen file from Synapse and collect unique drug names. + """ + query = ( + "select id, individualID, specimenID " + f"from {DRUG_SCREEN_TABLE} " + "where dataType='drug screen'" + ) + + files = syn.tableQuery(query).asDataFrame() + + logging.info("Reading %d drug-screen files for drug names", len(files)) + + names: set[str] = set() + + for _, row in files.iterrows(): + try: + entity = syn.get(row["id"]) + df = pd.read_csv(entity.path) + + except Exception as exc: + logging.warning("Skipping file %s: %s", row["id"], exc) + continue + + if "Drug" not in df.columns: + logging.warning("File %s has no 'Drug' column; skipping", row["id"]) + continue + + names.update( + df["Drug"] + .dropna() + .astype(str) + .str.strip() + .tolist() + ) + + names = sorted(n for n in names if n and n.lower() != "nan") + + logging.info("Collected %d unique drug names", len(names)) + + return names + + +# ---------------------------------------------------------------------------- +# Descriptor utility +# ---------------------------------------------------------------------------- + +def call_descriptor_table( + drug_file: str, + out_desc: str, + retries: int, + timeout: int, +) -> bool: + """ + Run the descriptor table builder as a subprocess. + + This function assumes build_drug_descriptor_table.py has a CLI accepting + --input and --output. + """ + cmd = [ + "python", + DESCRIPTOR_SCRIPT, + "--drugtable", + drug_file, + "--desctable", + out_desc, + ] + + ok = run_with_retries( + cmd, + retries, + timeout, + "descriptor_table", + ) + + if ok and not os.path.exists(out_desc): + logging.error( + "descriptor_table exited successfully but did not create %s", + out_desc, + ) + return False + + return ok + + +# ---------------------------------------------------------------------------- +# Main +# ---------------------------------------------------------------------------- + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + + parser.add_argument( + "--prev_drugs", + default="", + help="Comma-delimited list of existing drug TSV files for ID reuse.", + ) + + parser.add_argument( + "--out_drugs", + default=DEFAULT_OUT_DRUGS, + help=f"Output drugs TSV. Default: {DEFAULT_OUT_DRUGS}", + ) + + parser.add_argument( + "--out_desc", + default=DEFAULT_OUT_DESC, + help=f"Output drug descriptors TSV. Default: {DEFAULT_OUT_DESC}", + ) + + parser.add_argument( + "--skip-descriptors", + action="store_true", + help=( + "Skip the descriptor step entirely. Use when the descriptor " + "utility or its biomaRt/Ensembl dependency is failing. The drugs " + "table is still produced." + ), + ) + + parser.add_argument( + "--only-descriptors", + action="store_true", + help=( + "Only run the descriptor step, expecting --out_drugs to already " + "exist. Useful for retrying after a previous descriptor failure." + ), + ) + + parser.add_argument( + "--retries", + type=int, + default=DEFAULT_RETRIES, + help=f"Number of retry attempts. Default: {DEFAULT_RETRIES}", + ) + + parser.add_argument( + "--timeout", + type=int, + default=DEFAULT_TIMEOUT, + help=f"Timeout in seconds per attempt. Default: {DEFAULT_TIMEOUT}", + ) + + args = parser.parse_args() + + configure_logging() + + if args.skip_descriptors and args.only_descriptors: + logging.error("Cannot combine --skip-descriptors with --only-descriptors") + return 2 + + for path in (args.out_drugs, args.out_desc): + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + + # ------------------------------------------------------------------------ + # Drugs step + # ------------------------------------------------------------------------ + + if not args.only_descriptors: + syn = synapseclient.Synapse() + syn.login() + + names = collect_drug_names(syn) + + if not names: + logging.error("No drug names found; aborting") + return 1 + + prev = [ + p.strip() + for p in args.prev_drugs.split(",") + if p.strip() + ] + + ok = call_pubchem_retrieval( + names=names, + prev_drug_files=prev, + out_drugs=args.out_drugs, + retries=args.retries, + timeout=args.timeout, + ) + + if not ok: + logging.error("pubchem_retrieval failed; aborting") + return 1 + + if not os.path.exists(args.out_drugs): + logging.error( + "pubchem_retrieval did not produce %s", + args.out_drugs, + ) + return 1 + + logging.info("Wrote %s", args.out_drugs) + + # ------------------------------------------------------------------------ + # Descriptor step + # ------------------------------------------------------------------------ + + if args.skip_descriptors: + logging.warning( + "Skipping descriptor step. Re-run with --only-descriptors to " + "produce %s when descriptor dependencies are reachable.", + args.out_desc, + ) + return 0 + + if not os.path.exists(args.out_drugs): + logging.error( + "Drugs file %s missing; cannot run descriptors", + args.out_drugs, + ) + return 1 + + ok = call_descriptor_table( + drug_file=args.out_drugs, + out_desc=args.out_desc, + retries=args.retries, + timeout=args.timeout, + ) + + if not ok: + logging.error( + "Descriptor step failed after %d attempts. The drugs file %s is " + "complete; re-run with --only-descriptors when descriptor " + "dependencies are reachable to produce %s.", + args.retries, + args.out_drugs, + args.out_desc, + ) + + # Exit 0 because the primary drugs table was produced. Descriptors can + # be retried later. Use --skip-descriptors to suppress this explicitly. + return 0 + + logging.info("Wrote %s", args.out_desc) + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/coderbuild/cnf/04-experiments-cnf.py b/coderbuild/cnf/04-experiments-cnf.py new file mode 100644 index 00000000..fe1e216b --- /dev/null +++ b/coderbuild/cnf/04-experiments-cnf.py @@ -0,0 +1,478 @@ +#!/usr/bin/env python3 +""" +04-experiments-cnf.py: generate cnf_experiments.tsv. + +Discovers every drug-screen viability file under three Synapse parent folders, +attaches a specimen ID by reading the file's `specimenID` annotation (with a +filename-pattern fallback), then per (specimen, drug): + + * Multi-concentration measurements are run through the standard + coderbuild/utils/fit_curve.py utility, producing fit_auc / fit_ic50 / + fit_einf / fit_hs / fit_r2. + * Single-concentration measurements at 1 μM are recorded with + metric = 'uM_viability' and value = Viability_percentage / 100. + Schema dependency: 'uM_viability' must be in the ResponseMetric enum. + See schema_patch.md. + +Specimen attribution policy: + 1. Trust the file's `specimenID` annotation if present. + 2. Fall back to parsing the filename (e.g. "NF0017_T2_Viabilities.csv" + → "NF0017_T2"). + 3. Run the result through cnf_utils.classify_specimen as a safety net + (handles dot-separators, suffixes, case variants). Files whose + specimen can't be classified are skipped with a warning. + +Dual-mapping: one viability file's results should be attributed to two +specimens. Handled via SPECIMEN_DUAL_MAPPINGS at the top of the script; +add new entries as more cases come up. + +Usage: + python 04-experiments-cnf.py + [--output /tmp/cnf_experiments.tsv] + [--parents syn51301414,syn51301420,syn51301426] +""" + +import argparse +import logging +import os +import re +import subprocess +import sys +import tempfile + +import pandas as pd +import synapseclient + +from cnf_utils import classify_specimen + + +# ---------------------------------------------------------------------------- +# Constants +# ---------------------------------------------------------------------------- +# Synapse folders to walk for drug-screen viability files +DEFAULT_DRUG_SCREEN_PARENTS = [ + "syn51301414", + "syn51301420", + "syn51301426", +] + +TIME = 120 +TIME_UNIT = "hours" +STUDY = "cnf" +SOURCE = "Synapse" +SINGLE_DOSE_UM = 1.0 +SINGLE_DOSE_METRIC = "uM_viability" + +FIT_CURVE_SCRIPT = os.environ.get( + "CODERDATA_FIT_CURVE_SCRIPT", "/coderbuild/utils/fit_curve.py") + +DEFAULT_OUTPUT = "/tmp/cnf_experiments.tsv" + +# Regex to recover specimen from filenames like "NF0017_T2_Viabilities.csv" +# or "NF0021_T1_Onalespid_1uM_Viabilities.csv" +FILENAME_SPECIMEN_RE = re.compile( + r"^(NF\d{4}(?:_T\d+)?(?:_[A-Za-z0-9._-]+?)?)" # specimen prefix + r"_(?:Viabilities|viabilities|Viability|drug_screen)" # known label + r"\.(?:csv|tsv|txt)$", + re.IGNORECASE, +) + +# ---------------------------------------------------------------------------- +# DUAL-MAPPING TABLE +# ---------------------------------------------------------------------------- +# Maps a canonical specimen to any additional specimens that should receive +# the same drug-response rows. Add new entries as more cases come up. +# +# Known case: NF0021_T1_Onalespid_1uM was screened against the same drug +# panel as NF0021_T1 and should be available under both labels. +SPECIMEN_DUAL_MAPPINGS: dict[str, list[str]] = { + "NF0021_T1_Onalespid_1uM": ["NF0021_T1"], +} + + +def configure_logging(): + logging.basicConfig( + format="[%(asctime)s] [%(levelname)s] %(message)s", + datefmt="%Y-%m-%d %H:%M:%S", + level=logging.INFO, + ) + + +# ---------------------------------------------------------------------------- +# Specimen attribution +# ---------------------------------------------------------------------------- +def specimen_from_annotation(syn, file_id: str) -> str | None: + """Return the file's specimenID annotation, or None if missing.""" + try: + anno = syn.getAnnotations(file_id) + except Exception as exc: + logging.warning("getAnnotations failed for %s: %s", file_id, exc) + return None + val = anno.get("specimenID") + if isinstance(val, list): + val = val[0] if val else None + if val is None: + return None + s = str(val).strip() + return s or None + + +def specimen_from_filename(name: str) -> str | None: + """Parse 'NF0017_T2_Viabilities.csv' → 'NF0017_T2'.""" + m = FILENAME_SPECIMEN_RE.match(name) + if m: + return m.group(1) + return None + + +def resolve_specimen(syn, file_id: str, file_name: str) -> str | None: + """Resolve a file to a canonical specimen ID, trying annotation first. + + Returns canonical_id (e.g. "NF0017_T2") or None if both annotation and + filename parsing fail or the result doesn't classify as a known sample + type. + """ + raw = specimen_from_annotation(syn, file_id) + if not raw: + raw = specimen_from_filename(file_name) + if not raw: + return None + + # Safety-net classification — handles dot-separators, suffixes, case + result = classify_specimen(raw) + if result is None: + # Annotation or filename produced something that doesn't fit any + # known sample pattern. + logging.warning( + "Could not classify specimen '%s' (file %s, %s); skipping", + raw, file_id, file_name, + ) + return None + return result[0] + + +def expand_dual_mappings(specimen: str) -> list[str]: + """Return [specimen] + any additional specimens it's dual-mapped to.""" + return [specimen] + SPECIMEN_DUAL_MAPPINGS.get(specimen, []) + + +# ---------------------------------------------------------------------------- +# Synapse walking +# ---------------------------------------------------------------------------- +def walk_for_viability_files(syn, parent_id: str) -> list[dict]: + """Recursively yield File entities under a Synapse folder. + + Returns a list of {id, name} dicts for files whose name suggests they + contain viability data (drug-screen CSV/TSV). + """ + try: + import synapseutils + except ImportError: + logging.error("synapseutils not available; cannot walk %s", parent_id) + return [] + + files = [] + for dirpath, _subfolders, file_list in synapseutils.walk(syn, parent_id): + for fname, fid in file_list: + # Filter to plausible viability files. Accept anything ending in + # csv/tsv that has 'viab' or 'drug' in the name; fall back to + # 'csv' if neither matches but the file is small enough. + lower = fname.lower() + if not lower.endswith((".csv", ".tsv", ".txt")): + continue + if any(tok in lower for tok in ("viab", "drug", "screen")): + files.append({"id": fid, "name": fname, + "parent_path": dirpath}) + # Don't auto-grab arbitrary CSVs — we only want drug screens. + return files + + +def discover_drug_files(syn, parents: list[str]) -> pd.DataFrame: + """Walk all parent folders and dedupe by Synapse ID.""" + all_files = [] + for parent in parents: + logging.info("Walking %s ...", parent) + files = walk_for_viability_files(syn, parent) + logging.info(" found %d candidate viability files", len(files)) + all_files.extend(files) + + if not all_files: + return pd.DataFrame(columns=["id", "name"]) + + df = pd.DataFrame(all_files).drop_duplicates(subset="id").reset_index(drop=True) + logging.info("Discovered %d unique viability files across %d parents", + len(df), len(parents)) + return df + + +# ---------------------------------------------------------------------------- +# Read and combine drug-screen files +# ---------------------------------------------------------------------------- +def pull_screen_files(syn, parents: list[str]) -> pd.DataFrame: + """Walk parents, read each viability CSV, attach canonical specimen. + + Files whose specimen can't be resolved are skipped with a warning. + Files attributed to a dual-mapped specimen produce duplicate rows for + every additional specimen in the mapping. + """ + files_df = discover_drug_files(syn, parents) + if files_df.empty: + raise RuntimeError("No drug-screen files found under any parent") + + frames = [] + for _, row in files_df.iterrows(): + fid, fname = row["id"], row["name"] + + canonical = resolve_specimen(syn, fid, fname) + if not canonical: + continue + + try: + df = pd.read_csv(syn.get(fid).path) + except Exception as exc: + logging.warning("Skipping file %s: %s", fid, exc) + continue + + # Verify required columns + required = {"Drug", "Concentration_uM", "Viability_percentage"} + missing = required - set(df.columns) + if missing: + logging.warning("File %s missing %s; skipping", + fname, sorted(missing)) + continue + + # Attribute to canonical + any dual-mapped specimens + for specimen in expand_dual_mappings(canonical): + attributed = df.copy() + attributed["specimen"] = specimen + attributed["source_file"] = fid + frames.append(attributed) + + if not frames: + raise RuntimeError("No drug-screen files could be read successfully") + + combined = pd.concat(frames, ignore_index=True) + n_specimens = combined["specimen"].nunique() + n_drugs = combined["Drug"].nunique() + logging.info( + "Combined %d drug-screen rows across %d specimens × %d drugs", + len(combined), n_specimens, n_drugs, + ) + return combined + + +# ---------------------------------------------------------------------------- +# I/O helpers +# ---------------------------------------------------------------------------- +def load_sample_map(samples_path: str) -> dict[str, int]: + df = pd.read_csv(samples_path) + df["improve_sample_id"] = df["improve_sample_id"].astype(int) + return dict(zip(df["other_id"].astype(str), df["improve_sample_id"])) + + +def load_drug_map(drugs_path: str) -> dict[str, str]: + df = pd.read_csv(drugs_path, sep="\t") + if "chem_name" not in df.columns or "improve_drug_id" not in df.columns: + raise KeyError( + "drugs file missing required columns " + f"(have: {list(df.columns)})" + ) + mapping = { + str(name).strip().lower(): str(did) + for name, did in zip(df["chem_name"], df["improve_drug_id"]) + } + logging.info("Loaded %d drug name → improve_drug_id mappings", len(mapping)) + return mapping + + +def lookup_drug_id(name, drug_map: dict[str, str]) -> str | None: + if pd.isna(name): + return None + return drug_map.get(str(name).strip().lower()) + + +# ---------------------------------------------------------------------------- +# Single vs multi dose split +# ---------------------------------------------------------------------------- +def split_single_vs_multi(combined: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: + """Per (specimen, drug), separate single-dose vs multi-dose.""" + counts = ( + combined.groupby(["specimen", "Drug"])["Concentration_uM"] + .nunique() + .reset_index(name="n_conc") + ) + multi_keys = counts[counts["n_conc"] > 1] + single_keys = counts[counts["n_conc"] == 1] + + multi = combined.merge(multi_keys[["specimen", "Drug"]], + on=["specimen", "Drug"], how="inner") + single = combined.merge(single_keys[["specimen", "Drug"]], + on=["specimen", "Drug"], how="inner") + # Single-dose rows must be at the canonical 1 μM concentration. + single = single[single["Concentration_uM"] == SINGLE_DOSE_UM].copy() + + logging.info("Multi-dose rows: %d (across %d sample-drug pairs)", + len(multi), len(multi_keys)) + logging.info("Single-dose rows: %d (across %d sample-drug pairs at %g μM)", + len(single), len(single_keys), SINGLE_DOSE_UM) + return multi, single + + +# ---------------------------------------------------------------------------- +# Curve fitting (multi-dose) and single-dose formatting +# ---------------------------------------------------------------------------- +def run_curve_fit(multi: pd.DataFrame, sample_map, drug_map) -> pd.DataFrame: + if multi.empty: + return pd.DataFrame() + + multi = multi.copy() + multi["improve_sample_id"] = multi["specimen"].map(sample_map) + multi["improve_drug_id"] = multi["Drug"].apply( + lambda n: lookup_drug_id(n, drug_map)) + + n_drop = (multi["improve_sample_id"].isna() | + multi["improve_drug_id"].isna()).sum() + if n_drop: + logging.warning("Dropping %d multi-dose rows with unmapped sample/drug", + n_drop) + multi = multi.dropna(subset=["improve_sample_id", "improve_drug_id"]) + if multi.empty: + return pd.DataFrame() + + # fit_curve.py expects DOSE in μM and GROWTH as viability percentage + multi["DOSE"] = multi["Concentration_uM"].astype(float) + 1e-4 + multi["GROWTH"] = multi["Viability_percentage"].astype(float) + multi["time"] = TIME + multi["time_unit"] = TIME_UNIT + multi["study"] = STUDY + multi["source"] = SOURCE + + cols = ["DOSE", "GROWTH", "study", "source", "improve_sample_id", + "Drug", "time", "time_unit"] + multi_for_fit = multi[cols].copy() + multi_for_fit["improve_sample_id"] = ( + multi_for_fit["improve_sample_id"].astype(int)) + + with tempfile.TemporaryDirectory() as tmpdir: + in_path = os.path.join(tmpdir, "drug_response.tsv") + out_prefix = os.path.join(tmpdir, "cnf_curve") + multi_for_fit.to_csv(in_path, sep="\t", index=False) + + cmd = ["python", FIT_CURVE_SCRIPT, + "--input", in_path, "--output", out_prefix] + logging.info("Running fit_curve: %s", " ".join(cmd)) + subprocess.run(cmd, check=True) + + # fit_curve.py writes .0 in the reference snippet's behavior + out_path = out_prefix + ".0" + if not os.path.exists(out_path): + alt = out_prefix + if os.path.exists(alt): + out_path = alt + else: + raise FileNotFoundError( + f"fit_curve.py did not produce output at {out_path}" + ) + fitted = pd.read_csv(out_path, sep="\t") + + # fit_curve.py renames Drug → improve_drug_id, so the column is present + # but holds the raw drug name, not the SMI_ ID. Always remap through drug_map. + if "improve_drug_id" in fitted.columns: + name_col = "improve_drug_id" + elif "Drug" in fitted.columns: + name_col = "Drug" + else: + raise KeyError( + f"fit_curve output has neither improve_drug_id nor Drug " + f"column (have: {list(fitted.columns)})" + ) + fitted["improve_drug_id"] = fitted[name_col].apply( + lambda n: lookup_drug_id(n, drug_map)) + + fitted = fitted.dropna(subset=["improve_drug_id"]) + n_metrics = (fitted["dose_response_metric"].nunique() + if "dose_response_metric" in fitted.columns else 0) + logging.info("Curve-fit produced %d rows across %d metrics", + len(fitted), n_metrics) + return fitted + + +def format_single_dose(single: pd.DataFrame, sample_map, drug_map) -> pd.DataFrame: + if single.empty: + return pd.DataFrame() + + df = single.copy() + df["improve_sample_id"] = df["specimen"].map(sample_map) + df["improve_drug_id"] = df["Drug"].apply( + lambda n: lookup_drug_id(n, drug_map)) + + n_drop = (df["improve_sample_id"].isna() | + df["improve_drug_id"].isna()).sum() + if n_drop: + logging.warning("Dropping %d single-dose rows with unmapped sample/drug", + n_drop) + df = df.dropna(subset=["improve_sample_id", "improve_drug_id"]).copy() + + df["dose_response_metric"] = SINGLE_DOSE_METRIC + df["dose_response_value"] = df["Viability_percentage"].astype(float) / 100.0 + df["time"] = TIME + df["time_unit"] = TIME_UNIT + df["study"] = STUDY + df["source"] = SOURCE + df["improve_sample_id"] = df["improve_sample_id"].astype(int) + + cols = ["source", "improve_sample_id", "improve_drug_id", "study", + "time", "time_unit", "dose_response_metric", "dose_response_value"] + df = df[cols] + logging.info("Formatted %d single-dose rows", len(df)) + return df + + +# ---------------------------------------------------------------------------- +# Main +# ---------------------------------------------------------------------------- +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("samples", help="cnf_samples.csv from build_samples") + parser.add_argument("drugs", help="cnf_drugs.tsv from build_drugs") + parser.add_argument("--output", default=DEFAULT_OUTPUT) + parser.add_argument( + "--parents", default=",".join(DEFAULT_DRUG_SCREEN_PARENTS), + help="Comma-delimited Synapse folder IDs to walk for viability files", + ) + args = parser.parse_args() + + configure_logging() + + parents = [p.strip() for p in args.parents.split(",") if p.strip()] + if not parents: + logging.error("No parent folders provided") + sys.exit(2) + + sample_map = load_sample_map(args.samples) + drug_map = load_drug_map(args.drugs) + + syn = synapseclient.Synapse() + syn.login() + + combined = pull_screen_files(syn, parents) + multi, single = split_single_vs_multi(combined) + + fitted = run_curve_fit(multi, sample_map, drug_map) + single_df = format_single_dose(single, sample_map, drug_map) + + final = pd.concat([fitted, single_df], ignore_index=True) + cols = ["source", "improve_sample_id", "improve_drug_id", "study", + "time", "time_unit", "dose_response_metric", "dose_response_value"] + for c in cols: + if c not in final.columns: + final[c] = pd.NA + final = final[cols] + + os.makedirs(os.path.dirname(args.output) or ".", exist_ok=True) + final.to_csv(args.output, sep="\t", index=False) + logging.info("Wrote %s (%d rows)", args.output, len(final)) + + +if __name__ == "__main__": + sys.exit(main() or 0) diff --git a/coderbuild/cnf/build_drugs.sh b/coderbuild/cnf/build_drugs.sh new file mode 100644 index 00000000..9402f178 --- /dev/null +++ b/coderbuild/cnf/build_drugs.sh @@ -0,0 +1,20 @@ +#!/usr/bin/env bash +# build_drugs.sh - wraps 03-drugs-cnf.py per coderdata convention. +# +# Usage: +# build_drugs.sh +# +# is a comma-delimited list of drug TSV files from +# previous datasets in the build sequence. Their improve_drug_id values +# get reused when the canonical SMILES match. + +set -euo pipefail + +PREV_DRUGS="${1:-}" + +python 03-drugs-cnf.py \ + --prev_drugs "$PREV_DRUGS" \ + --out_drugs /tmp/cnf_drugs.tsv \ + --out_desc /tmp/cnf_drug_descriptors.tsv + +python build_drug_desc.py --drugtable /tmp/cnf_drugs.tsv --desctable /tmp/cnf_drug_descriptors.tsv.gz diff --git a/coderbuild/cnf/build_exp.sh b/coderbuild/cnf/build_exp.sh new file mode 100644 index 00000000..18c8b107 --- /dev/null +++ b/coderbuild/cnf/build_exp.sh @@ -0,0 +1,12 @@ +#!/usr/bin/env bash +# build_exp.sh - wraps 04-experiments-cnf.py per coderdata convention. +# +# Usage: +# build_exp.sh + +set -euo pipefail + +SAMPLES="${1:?Usage: build_exp.sh }" +DRUGS="${2:?Usage: build_exp.sh }" + +python 04-experiments-cnf.py "$SAMPLES" "$DRUGS" --output /tmp/cnf_experiments.tsv \ No newline at end of file