diff --git a/DESCRIPTION b/DESCRIPTION index d86ec6a..e3e40c1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -31,7 +31,8 @@ Imports: dplyr, stats, parallel, - data.table + data.table, + plotly Suggests: MSstats, MSstatsTMT, diff --git a/NAMESPACE b/NAMESPACE index d7a1451..24d4bd3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,8 +4,10 @@ export(MSstatsPrepareDoseResponseFit) export(convertGroupToNumericDose) export(doseResponseFit) export(futureExperimentSimulation) +export(plot_tpr_power_curve) export(predictIC50) export(predictIC50Parallel) +export(run_tpr_simulation) export(visualizeResponseProtein) import(dplyr) importFrom(BiocParallel,bplapply) @@ -15,6 +17,7 @@ importFrom(dplyr,arrange) importFrom(dplyr,distinct) importFrom(dplyr,filter) importFrom(dplyr,group_by) +importFrom(dplyr,if_else) importFrom(dplyr,mutate) importFrom(dplyr,select) importFrom(dplyr,summarise) @@ -35,6 +38,7 @@ importFrom(ggplot2,scale_fill_manual) importFrom(ggplot2,scale_x_continuous) importFrom(ggplot2,scale_y_continuous) importFrom(ggplot2,theme) +importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_classic) importFrom(ggplot2,theme_minimal) importFrom(ggplot2,ylim) @@ -42,6 +46,8 @@ importFrom(parallel,clusterExport) importFrom(parallel,makeCluster) importFrom(parallel,parLapply) importFrom(parallel,stopCluster) +importFrom(plotly,ggplotly) +importFrom(plotly,layout) importFrom(stats,approx) importFrom(stats,p.adjust) importFrom(stats,pf) diff --git a/R/ExperimentalDesignSimulation.R b/R/ExperimentalDesignSimulation.R index e05f794..41bfc1d 100644 --- a/R/ExperimentalDesignSimulation.R +++ b/R/ExperimentalDesignSimulation.R @@ -222,23 +222,30 @@ futureExperimentSimulation = function(N_proteins = 300, )) } - # Get mean profile across proteins + # Get mean profile across proteins using string-based grouping to avoid floating point precision issues profile = drug_data %>% filter(protein %in% protein_ids) %>% - mutate(dose_nM = round(dose * 1e9)) %>% # Convert M to nM and round to handle precision - group_by(dose_nM) %>% + mutate(dose_str = as.character(dose)) %>% + group_by(dose_str) %>% summarise( LogIntensities = mean(response, na.rm = TRUE), + dose_numeric = dose[1], .groups = 'drop' - ) %>% - rename(dose = dose_nM) + ) - # Create complete profile with all concentrations + # Match profile doses to target concentrations using nearest-neighbor on log scale. + # This handles unit mismatches without assuming a fixed conversion factor. complete_profile = data.frame(dose = concentrations) - # Merge with existing data + if (nrow(profile) > 0 && !any(profile$dose_numeric != 0)) { + stop("No non-zero dose measurements found for the selected protein.") + } + + # Match by string representation to avoid floating point issues complete_profile = complete_profile %>% - left_join(profile, by = "dose") + mutate(dose_str = as.character(dose)) %>% + left_join(profile %>% select(dose_str, LogIntensities), by = "dose_str") %>% + select(-dose_str) # Fill missing values with interpolation or nearest neighbor for (i in which(is.na(complete_profile$LogIntensities))) { @@ -306,7 +313,7 @@ futureExperimentSimulation = function(N_proteins = 300, return(complete_profile) } - # Extract templates for each category + # Extract templates: use user data for specified proteins, baseline for others template = list( strong_interaction = .getMeanProfile(strong_proteins, drug_data, concentrations), weak_interaction = .getMeanProfile(weak_proteins, drug_data, concentrations), diff --git a/R/TPR_Power_Curve.R b/R/TPR_Power_Curve.R new file mode 100644 index 0000000..8f066e6 --- /dev/null +++ b/R/TPR_Power_Curve.R @@ -0,0 +1,314 @@ +# ============================================================================ +# TPR Power Curve - Simulation and Visualization +# ============================================================================ + +#' Build concentration subsets for dose-response power analysis +#' +#' Given a user's full set of experimental concentrations, generates subsets +#' of size 2 through N by greedily selecting doses that are farthest from +#' the log-median of the already-selected set. Control (0) and the highest +#' dose are always included. +#' +#' @param concentrations Numeric vector of all available concentrations. +#' @param dose_range Integer vector of length 2, c(min, max) specifying the +#' range of subset sizes to generate. +#' +#' @return A named list where each element is a numeric vector of concentrations. +#' Names are the subset sizes (e.g., "2", "3", ...). +#' @noRd +.build_concentration_ladders <- function(concentrations, dose_range) { + concentrations <- sort(unique(concentrations)) + if (length(concentrations) < 2) { + stop("At least 2 unique concentrations are required.") + } + if (!(0 %in% concentrations)) { + stop("Concentrations must include 0 (control).") + } + + if (dose_range[2] > length(concentrations)) { + stop(paste0("dose_range maximum (", dose_range[2], + ") exceeds the number of available concentrations (", + length(concentrations), ").")) + } + + control <- min(concentrations) + highest <- max(concentrations) + candidates <- setdiff(concentrations, c(control, highest)) + + # Greedily build up from {control, highest} by picking the dose + # farthest from log(median) of the current set at each step + selected <- c(control, highest) + selection_order <- list() + selection_order[["2"]] <- sort(selected) + + remaining <- candidates + max_k <- min(dose_range[2], length(concentrations)) + + for (k in 3:max_k) { + if (length(remaining) == 0) break + + log_selected <- log10(selected + 1) + log_median <- median(log_selected) + log_remaining <- log10(remaining + 1) + + distances <- abs(log_remaining - log_median) + best_idx <- which.max(distances) + pick <- remaining[best_idx] + + selected <- c(selected, pick) + remaining <- remaining[-best_idx] + selection_order[[as.character(k)]] <- sort(selected) + } + + # Filter to requested dose_range + k_range <- seq(dose_range[1], dose_range[2]) + result <- selection_order[as.character(k_range)] + return(result) +} + +#' Run a single TPR simulation for one replicate/concentration combination +#' +#' @param n_rep Integer. Number of replicates per dose. +#' @param concs Numeric vector. Concentrations to use. +#' @param n_proteins Integer. Number of proteins to simulate. +#' @param data Prepared dose-response data for user-based templates. +#' @param protein Protein ID for strong interaction template. +#' @param seed Integer. Base seed for reproducibility. +#' +#' @return A data.frame with columns: Interaction, TPR, N_rep, NumConcs. +#' @noRd +.run_one_tpr_simulation <- function(n_rep, concs, n_proteins, data, + protein, seed = 123) { + set.seed(seed) + + sim_args <- list( + N_proteins = n_proteins, + N_rep = n_rep, + Concentrations = concs, + IC50_Prediction = FALSE + ) + sim_args$data <- data + sim_args$strong_proteins <- protein + + temp_res <- do.call(futureExperimentSimulation, sim_args) + temp_res$Hit_Rates_Data |> + dplyr::filter(Category == "TPR (Strong)") |> + dplyr::mutate( + N_rep = n_rep, + NumConcs = length(concs), + Interaction = dplyr::if_else(Category == "TPR (Strong)", "Strong", "Weak") + ) |> + dplyr::select(Interaction, TPR = Percent, N_rep, NumConcs) +} + +#' Simulate detection power across experimental design configurations +#' +#' Evaluates how well a dose-response proteomics experiment can detect +#' drug-protein interactions under different experimental designs. For each +#' combination of replicate count and number of doses, the function simulates +#' a full experiment and measures the true positive rate (TPR); the +#' percentage of known interacting proteins that are correctly identified. +#' +#' This helps experimentalists answer practical questions like: "If I can +#' only afford 3 replicates per dose, how many dose levels do I need to +#' reliably detect weak interactions?" +#' +#' Concentration subsets are selected automatically from the user's available +#' doses: control (0) and the highest dose are always included, and +#' intermediate doses are added by choosing the one farthest from the +#' log-median of the already-selected set. +#' +#' @param rep_range Integer vector of length 2, \code{c(min, max)}. The range +#' of biological replicates per dose to evaluate. Each value in the sequence +#' \code{min:max} produces one line in the resulting power curve. +#' @param concentrations Numeric vector. The full set of drug concentrations +#' available in the experiment (e.g., \code{c(0, 1, 3, 10, 30, 100, 300, 1000, 3000)}). +#' @param dose_range Integer vector of length 2, \code{c(min, max)}. The range +#' of dose counts to sweep. For example, \code{c(2, 9)} evaluates designs +#' with 2 doses, 3 doses, ..., up to 9 doses. +#' @param data User's prepared dose-response data (e.g., from +#' \code{MSstatsPrepareDoseResponseFit}). Protein-specific templates are +#' extracted from this data for the simulation. +#' @param protein Character string specifying a protein ID to use as +#' the interaction template for the power analysis. +#' @param n_proteins Integer. Number of proteins to simulate per run. Larger +#' values give more stable estimates but increase runtime. Default: 1000. +#' +#' @return A \code{data.frame} with columns: +#' \describe{ +#' \item{Interaction}{Character. "Strong" or "Weak".} +#' \item{TPR}{Numeric. True positive rate as a percentage (0-100).} +#' \item{N_rep}{Integer. Number of replicates used in this simulation.} +#' \item{NumConcs}{Integer. Number of concentrations used in this simulation.} +#' } +#' +#' @examples +#' \dontrun{ +#' # Prepare user data +#' mock_data <- data.frame( +#' protein = rep("P1", 6), +#' drug = c("DMSO", "DMSO", "Drug1", "Drug1", "Drug1", "Drug1"), +#' dose = c(0, 0, 10, 100, 1000, 1000), +#' response = c(20, 20, 18, 15, 12, 12) +#' ) +#' +#' # Run power analysis +#' results <- run_tpr_simulation( +#' rep_range = c(1, 3), +#' concentrations = c(0, 10, 100, 1000), +#' dose_range = c(2, 4), +#' data = mock_data, +#' protein = "P1", +#' n_proteins = 300 +#' ) +#' +#' # Visualize results +#' plot_tpr_power_curve(results) +#' } +#' +#' @importFrom dplyr filter mutate select if_else +#' @export +run_tpr_simulation <- function(rep_range, concentrations, dose_range, + data, protein, n_proteins = 1000) { + if (!is.numeric(rep_range) || length(rep_range) != 2 || rep_range[1] > rep_range[2]) { + stop("rep_range must be a numeric vector of length 2 with c(min, max) where min <= max.") + } + if (!is.numeric(dose_range) || length(dose_range) != 2 || dose_range[1] > dose_range[2]) { + stop("dose_range must be a numeric vector of length 2 with c(min, max) where min <= max.") + } + if (dose_range[1] < 2) { + stop("dose_range minimum must be at least 2 (control + one treatment).") + } + + conc_subsets <- .build_concentration_ladders(concentrations, dose_range) + k_grid <- as.integer(names(conc_subsets)) + rep_grid <- seq(rep_range[1], rep_range[2]) + + grid_df <- expand.grid(N_rep = rep_grid, k_conc = k_grid) + results <- do.call(rbind, lapply(seq_len(nrow(grid_df)), function(i) { + k <- as.character(grid_df$k_conc[i]) + concs <- conc_subsets[[k]] + tryCatch( + .run_one_tpr_simulation( + n_rep = grid_df$N_rep[i], + concs = concs, + n_proteins = n_proteins, + data = data, + protein = protein + ), + error = function(e) { + warning(paste("Simulation failed for N_rep=", grid_df$N_rep[i], + ", k_conc=", grid_df$k_conc[i], ":", conditionMessage(e))) + NULL + } + ) + })) + + if (is.null(results) || nrow(results) == 0) { + stop("All simulation runs failed. Please check your input parameters.") + } + return(results) +} + +#' Create a single TPR panel plot with color gradient by replicate count +#' +#' @param data Data frame with TPR results for one interaction type. +#' @param k_grid Integer vector. X-axis breaks. +#' @param show_legend Logical. Whether to display the legend. +#' +#' @return A ggplot object. +#' @importFrom ggplot2 ggplot aes geom_line geom_point scale_x_continuous +#' scale_y_continuous scale_color_manual labs theme_bw theme element_text +#' @noRd +.make_tpr_panel <- function(data, k_grid, show_legend = FALSE) { + rep_levels <- sort(unique(data$N_rep)) + n_levels <- length(rep_levels) + + # 10 distinct colors for up to 10 replicate levels + all_colors <- c( + "#1f78b4", # strong blue + "#ff964f", # orange + "#17becf", # cyan + "#1b9e77", # teal + "#f0b6d5", # pink + "#DC143C", # crimson + "#6a3d9a", # purple + "#9e9ac8", # lavender + "#525252", # cool dark gray + "#08306b" # deep navy + ) + color_palette <- setNames(all_colors[seq_len(n_levels)], as.character(rep_levels)) + + p <- ggplot2::ggplot(data, + ggplot2::aes(x = NumConcs, y = TPR, + color = factor(N_rep), group = factor(N_rep))) + + ggplot2::geom_line(linewidth = 1.2) + + ggplot2::geom_point(size = 2) + + ggplot2::scale_x_continuous(breaks = k_grid) + + ggplot2::scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20)) + + ggplot2::scale_color_manual(values = color_palette) + + ggplot2::labs( + x = "Number of concentrations", + y = "True Positive Rate (%)", + color = "Replicates" + ) + + ggplot2::theme_bw(base_size = 14) + + ggplot2::theme( + plot.title = ggplot2::element_text(face = "bold", hjust = 0.5) + ) + + if (!show_legend) { + p <- p + ggplot2::theme(legend.position = "none") + } + p +} + +#' Visualize detection power across experimental designs +#' +#' Creates an interactive plot showing how the true positive rate +#' (detection power) changes with the number of doses and replicates. +#' Only the results for the user-selected protein template (passed as +#' the strong interaction category) are displayed. Line color shading goes +#' from light (fewest replicates) to dark (most replicates). +#' +#' @param static Logical. If \code{TRUE}, returns a static \code{ggplot} +#' object suitable for PDF export. Default: \code{FALSE}. +#' +#' @return An interactive \code{plotly} object, or a static \code{ggplot} +#' object if \code{static = TRUE}. +#' +#' @examples +#' \dontrun{ +#' results <- run_tpr_simulation( +#' rep_range = c(1, 3), +#' concentrations = c(0, 10, 100, 1000), +#' dose_range = c(2, 4), +#' n_proteins = 300 +#' ) +#' plot_tpr_power_curve(results) +#' } +#' +#' @importFrom ggplot2 ggplot aes geom_line geom_point scale_x_continuous +#' scale_y_continuous scale_color_manual labs theme_bw theme element_text +#' @importFrom plotly ggplotly layout +#' @export +plot_tpr_power_curve <- function(simulation_results, static = FALSE) { + k_grid <- sort(unique(simulation_results$NumConcs)) + + # Use only the Strong interaction results (user-selected protein template) + results_protein <- simulation_results[simulation_results$Interaction == "Strong", ] + + p <- .make_tpr_panel(results_protein, k_grid, TRUE) + + ggplot2::ggtitle("Interaction detection power") + + if (static) return(p) + + plotly::ggplotly(p) |> plotly::layout( + margin = list(t = 60), + annotations = list( + list(text = "Interaction detection power", x = 0.5, y = 1.08, + xref = "paper", yref = "paper", showarrow = FALSE, + font = list(size = 14), xanchor = "center") + ) + ) +} \ No newline at end of file diff --git a/man/DoseResponseFit.Rd b/man/DoseResponseFit.Rd index 07dd21f..3033264 100644 --- a/man/DoseResponseFit.Rd +++ b/man/DoseResponseFit.Rd @@ -66,7 +66,7 @@ interaction_results <- doseResponseFit( print(interaction_results) # Check significant interactions -significant <- interaction_results[interaction_results$adjust_pval < 0.05, ] +significant <- interaction_results[interaction_results$adj.pvalue < 0.05, ] print(paste("Found", nrow(significant), "significant interactions")) \dontrun{ diff --git a/man/plot_tpr_power_curve.Rd b/man/plot_tpr_power_curve.Rd new file mode 100644 index 0000000..b3c74d4 --- /dev/null +++ b/man/plot_tpr_power_curve.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/TPR_Power_Curve.R +\name{plot_tpr_power_curve} +\alias{plot_tpr_power_curve} +\title{Visualize detection power across experimental designs} +\usage{ +plot_tpr_power_curve(simulation_results, static = FALSE) +} +\arguments{ +\item{static}{Logical. If \code{TRUE}, returns a static \code{ggplot} +object suitable for PDF export. Default: \code{FALSE}.} +} +\value{ +An interactive \code{plotly} object, or a static \code{ggplot} +object if \code{static = TRUE}. +} +\description{ +Creates an interactive plot showing how the true positive rate +(detection power) changes with the number of doses and replicates. +Only the results for the user-selected protein template (passed as +the strong interaction category) are displayed. Line color shading goes +from light (fewest replicates) to dark (most replicates). +} +\examples{ +\dontrun{ +results <- run_tpr_simulation( + rep_range = c(1, 3), + concentrations = c(0, 10, 100, 1000), + dose_range = c(2, 4), + n_proteins = 300 +) +plot_tpr_power_curve(results) +} + +} diff --git a/man/run_tpr_simulation.Rd b/man/run_tpr_simulation.Rd new file mode 100644 index 0000000..db7a1c2 --- /dev/null +++ b/man/run_tpr_simulation.Rd @@ -0,0 +1,88 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/TPR_Power_Curve.R +\name{run_tpr_simulation} +\alias{run_tpr_simulation} +\title{Simulate detection power across experimental design configurations} +\usage{ +run_tpr_simulation( + rep_range, + concentrations, + dose_range, + data, + protein, + n_proteins = 1000 +) +} +\arguments{ +\item{rep_range}{Integer vector of length 2, \code{c(min, max)}. The range +of biological replicates per dose to evaluate. Each value in the sequence +\code{min:max} produces one line in the resulting power curve.} + +\item{concentrations}{Numeric vector. The full set of drug concentrations +available in the experiment (e.g., \code{c(0, 1, 3, 10, 30, 100, 300, 1000, 3000)}).} + +\item{dose_range}{Integer vector of length 2, \code{c(min, max)}. The range +of dose counts to sweep. For example, \code{c(2, 9)} evaluates designs +with 2 doses, 3 doses, ..., up to 9 doses.} + +\item{data}{User's prepared dose-response data (e.g., from +\code{MSstatsPrepareDoseResponseFit}). Protein-specific templates are +extracted from this data for the simulation.} + +\item{protein}{Character string specifying a protein ID to use as +the interaction template for the power analysis.} + +\item{n_proteins}{Integer. Number of proteins to simulate per run. Larger +values give more stable estimates but increase runtime. Default: 1000.} +} +\value{ +A \code{data.frame} with columns: +\describe{ +\item{Interaction}{Character. "Strong" or "Weak".} +\item{TPR}{Numeric. True positive rate as a percentage (0-100).} +\item{N_rep}{Integer. Number of replicates used in this simulation.} +\item{NumConcs}{Integer. Number of concentrations used in this simulation.} +} +} +\description{ +Evaluates how well a dose-response proteomics experiment can detect +drug-protein interactions under different experimental designs. For each +combination of replicate count and number of doses, the function simulates +a full experiment and measures the true positive rate (TPR); the +percentage of known interacting proteins that are correctly identified. +} +\details{ +This helps experimentalists answer practical questions like: "If I can +only afford 3 replicates per dose, how many dose levels do I need to +reliably detect weak interactions?" + +Concentration subsets are selected automatically from the user's available +doses: control (0) and the highest dose are always included, and +intermediate doses are added by choosing the one farthest from the +log-median of the already-selected set. +} +\examples{ +\dontrun{ +# Prepare user data +mock_data <- data.frame( + protein = rep("P1", 6), + drug = c("DMSO", "DMSO", "Drug1", "Drug1", "Drug1", "Drug1"), + dose = c(0, 0, 10, 100, 1000, 1000), + response = c(20, 20, 18, 15, 12, 12) +) + +# Run power analysis +results <- run_tpr_simulation( + rep_range = c(1, 3), + concentrations = c(0, 10, 100, 1000), + dose_range = c(2, 4), + data = mock_data, + protein = "P1", + n_proteins = 300 +) + +# Visualize results +plot_tpr_power_curve(results) +} + +} diff --git a/tests/testthat/test-ExperimentalDesignSimulation.R b/tests/testthat/test-ExperimentalDesignSimulation.R index 95963fc..8e85747 100644 --- a/tests/testthat/test-ExperimentalDesignSimulation.R +++ b/tests/testthat/test-ExperimentalDesignSimulation.R @@ -133,12 +133,12 @@ test_that(".extractTemplatesFromData handles missing doses via interpolation", { weak_proteins = NULL, no_interaction_proteins = NULL, drug_name = "Drug1", - concentrations = c(0, 10, 100, 1000, 10000) # More concentrations than data + concentrations = c(0, 10, 100, 1000, 3000) # All within default template range ) # Should have all requested concentrations expect_equal(nrow(template$strong_interaction), 5) - expect_equal(template$strong_interaction$dose, c(0, 10, 100, 1000, 10000)) + expect_equal(template$strong_interaction$dose, c(0, 10, 100, 1000, 3000)) }) test_that(".extractTemplatesFromData handles empty protein lists", { @@ -158,7 +158,12 @@ test_that(".extractTemplatesFromData handles empty protein lists", { concentrations = c(0, 1000) ) - # Null protein lists should get flat profiles + # Null protein lists get flat baseline profiles + expect_equal(nrow(template$weak_interaction), 2) + expect_equal(nrow(template$no_interaction), 2) + expect_true(all(c("dose", "Intensity", "LogIntensities", "ratio") %in% names(template$weak_interaction))) + expect_true(all(c("dose", "Intensity", "LogIntensities", "ratio") %in% names(template$no_interaction))) + # Flat profiles should have ratio = 1 for all doses expect_equal(template$weak_interaction$ratio, c(1, 1)) expect_equal(template$no_interaction$ratio, c(1, 1)) }) diff --git a/tests/testthat/test-TPR_Power_Curve.R b/tests/testthat/test-TPR_Power_Curve.R new file mode 100644 index 0000000..35c0b9b --- /dev/null +++ b/tests/testthat/test-TPR_Power_Curve.R @@ -0,0 +1,146 @@ +test_that("build_concentration_ladders returns correct subset sizes", { + concs <- c(0, 1, 3, 10, 30, 100, 300, 1000, 3000) + result <- MSstatsResponse:::.build_concentration_ladders(concs, c(2, 5)) + + expect_equal(length(result), 4) + expect_equal(names(result), c("2", "3", "4", "5")) + + for (k in names(result)) { + expect_true(0 %in% result[[k]]) + expect_true(3000 %in% result[[k]]) + expect_equal(length(result[[k]]), as.integer(k)) + } +}) + +test_that("build_concentration_ladders always includes control and max", { + concs <- c(0, 5, 50, 200, 1000) + result <- MSstatsResponse:::.build_concentration_ladders(concs, c(2, 5)) + + for (k in names(result)) { + expect_true(0 %in% result[[k]]) + expect_true(1000 %in% result[[k]]) + } +}) + +test_that("build_concentration_ladders errors on fewer than 2 concentrations", { + expect_error( + MSstatsResponse:::.build_concentration_ladders(c(0), c(2, 3)), + "At least 2 unique" + ) +}) + +test_that("build_concentration_ladders errors without control dose", { + expect_error( + MSstatsResponse:::.build_concentration_ladders(c(1, 10, 100), c(2, 3)), + "must include 0" + ) +}) + +test_that("build_concentration_ladders handles arbitrary user concentrations", { + concs <- c(0, 2, 4, 20, 50, 100, 500, 1000) + result <- MSstatsResponse:::.build_concentration_ladders(concs, c(2, 6)) + + expect_equal(length(result), 5) + for (k in names(result)) { + subset <- result[[k]] + expect_true(all(subset %in% concs)) + } +}) + +test_that("run_tpr_simulation validates rep_range", { + mock_data <- data.frame(protein = "P1", drug = "D1", dose = 0, response = 20) + expect_error( + run_tpr_simulation(c(3, 1), c(0, 100, 1000), c(2, 3), data = mock_data, protein = "P1"), + "min <= max" + ) + expect_error( + run_tpr_simulation("bad", c(0, 100, 1000), c(2, 3), data = mock_data, protein = "P1"), + "numeric vector" + ) +}) + +test_that("run_tpr_simulation validates dose_range", { + mock_data <- data.frame(protein = "P1", drug = "D1", dose = 0, response = 20) + expect_error( + run_tpr_simulation(c(1, 2), c(0, 100, 1000), c(5, 2), data = mock_data, protein = "P1"), + "min <= max" + ) + expect_error( + run_tpr_simulation(c(1, 2), c(0, 100, 1000), c(1, 3), data = mock_data, protein = "P1"), + "at least 2" + ) +}) + +test_that("run_tpr_simulation returns correct structure", { + mock_data <- data.frame( + protein = rep("P1", 6), + drug = c("DMSO", "DMSO", "Drug1", "Drug1", "Drug1", "Drug1"), + dose = c(0, 0, 10, 100, 1000, 1000), + response = c(20, 20, 18, 15, 12, 12) + ) + + results <- run_tpr_simulation( + rep_range = c(1, 3), + concentrations = c(0, 10, 100, 1000), + dose_range = c(2, 4), + data = mock_data, + protein = "P1", + n_proteins = 200 + ) + + expect_true(is.data.frame(results)) + expect_true(all(c("Interaction", "TPR", "N_rep", "NumConcs") %in% colnames(results))) + expect_true(all(results$Interaction == "Strong")) + expect_true(all(results$TPR >= 0 & results$TPR <= 100)) + + noise_margin <- 5 + + for (nc in unique(results$NumConcs)) { + subset <- results[results$NumConcs == nc, ] + subset <- subset[order(subset$N_rep), ] + if (nrow(subset) >= 2) { + for (j in 2:nrow(subset)) { + expect_true(subset$TPR[j] >= subset$TPR[j - 1] - noise_margin, + info = paste0("TPR should not drop significantly between rep ", + subset$N_rep[j - 1], " and rep ", subset$N_rep[j], + " at NumConcs = ", nc, + " (was ", subset$TPR[j - 1], ", got ", subset$TPR[j], ")")) + } + } + } + + for (nr in unique(results$N_rep)) { + subset <- results[results$N_rep == nr, ] + subset <- subset[order(subset$NumConcs), ] + if (nrow(subset) >= 2) { + for (j in 2:nrow(subset)) { + expect_true(subset$TPR[j] >= subset$TPR[j - 1] - noise_margin, + info = paste0("TPR should not drop significantly between ", + subset$NumConcs[j - 1], " and ", subset$NumConcs[j], + " concentrations at N_rep = ", nr, + " (was ", subset$TPR[j - 1], ", got ", subset$TPR[j], ")")) + } + } + } +}) + +test_that("plot_tpr_power_curve returns a plotly object", { + mock_data <- data.frame( + protein = rep("P1", 6), + drug = c("DMSO", "DMSO", "Drug1", "Drug1", "Drug1", "Drug1"), + dose = c(0, 0, 10, 100, 1000, 1000), + response = c(20, 20, 18, 15, 12, 12) + ) + + results <- run_tpr_simulation( + rep_range = c(1, 3), + concentrations = c(0, 10, 100, 1000), + dose_range = c(2, 4), + data = mock_data, + protein = "P1", + n_proteins = 200 + ) + + p <- plot_tpr_power_curve(results) + expect_true(inherits(p, "plotly")) +}) \ No newline at end of file