## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----install, eval=FALSE------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("MSstatsResponse")
# 

## ----load_libraries, message=FALSE--------------------------------------------
library(MSstatsResponse)
library(dplyr)
library(ggplot2)
library(data.table)

# Optional: for upstream data processing
# library(MSstats)
# library(MSstatsTMT)

## ----msstats_preprocessing, eval=FALSE----------------------------------------
# # Read raw data (example with Spectronaut output)
# raw_data <- read_tsv("path/to/spectronaut_report.tsv")
# 
# # Convert to MSstats format
# msstats_data <- SpectronauttoMSstatsFormat(raw_data)
# 
# # Process data: normalization and protein summarization
# processed_data <- dataProcess(
#   msstats_data,
#   normalization = "equalizeMedians",
#   summaryMethod = "TMP",
#   MBimpute = TRUE,
#   maxQuantileforCensored = 0.999
# )
# 
# # Extract protein-level data for dose-response analysis
# protein_level_data <- processed_data$ProteinLevelData

## ----load_example_data--------------------------------------------------------
data("DIA_MSstats_Normalized", package = "MSstatsResponse")
dia_normalized <- DIA_MSstats_Normalized

str(dia_normalized$ProteinLevelData[1:5, ])

## ----convert_doses------------------------------------------------------------
dose_info <- convertGroupToNumericDose(dia_normalized$ProteinLevelData$GROUP)

dia_normalized$ProteinLevelData <- dia_normalized$ProteinLevelData %>%
  mutate(
    dose_nM = dose_info$dose_nM,
    drug = dose_info$drug
  )

dia_normalized$ProteinLevelData %>%
  select(Protein, GROUP, drug, dose_nM) %>%
  head(10)

## ----prepare_data-------------------------------------------------------------
dia_prepared <- MSstatsPrepareDoseResponseFit(
  data = dia_normalized$ProteinLevelData,
  dose_column = "dose_nM",
  drug_column = "drug",
  protein_column = "Protein",
  log_abundance_column = "LogIntensities",
  transform_nM_to_M = TRUE
)

str(dia_prepared)

dia_prepared %>%
  filter(protein %in% c("PROTEIN_A", "PROTEIN_B")) %>%
  arrange(protein, drug, dose) %>%
  head(20)

## ----fit_dose_response, message=FALSE, warning=FALSE--------------------------
response_results <- doseResponseFit(
  data = dia_prepared,
  increasing = FALSE,
  transform_dose = TRUE,
  ratio_response = FALSE
)

response_results %>%
  select(Protein, drug, F_statistic, log2FC, adj.pvalue) %>%
  arrange(adj.pvalue) %>%
  head(10)

## ----summarize_results--------------------------------------------------------
n_total <- nrow(response_results)
n_significant <- sum(response_results$adj.pvalue < 0.05)

cat("Total protein-drug pairs tested:", n_total, "\n")
cat("Significant interactions (FDR < 0.05):", n_significant, "\n")
cat("Percentage significant:", round(100 * n_significant / n_total, 1), "%\n")

top_hits <- response_results %>%
  filter(adj.pvalue < 0.05) %>%
  arrange(adj.pvalue) %>%
  head(5)

print(top_hits)

## ----predict_ic50, message=FALSE, warning=FALSE-------------------------------
ic50_predictions <- predictIC50(
  data = dia_prepared,
  ratio_response = TRUE,
  transform_dose = TRUE,
  increasing = FALSE,
  bootstrap = TRUE,
  n_samples = 1000,
  alpha = 0.10
)

ic50_predictions %>%
  arrange(IC50) %>%
  head(10)

## ----parallel_ic50, eval=FALSE------------------------------------------------
# library(BiocParallel)
# if (.Platform$OS.type == "windows") {
#   register(SnowParam(workers = 4, type = "SOCK"))
# } else {
#   register(MulticoreParam(workers = 4))
# }
# 
# ic50_parallel <- predictIC50(
#   data = dia_prepared,
#   ratio_response = TRUE,
#   transform_dose = TRUE,
#   bootstrap = TRUE,
#   n_samples = 1000,
#   BPPARAM = bpparam()
# )

## ----visualize_single, fig.height=5, fig.width=7------------------------------
visualizeResponseProtein(
  data = dia_prepared,
  protein_name = "PROTEIN_A",
  drug_name = "Drug1",
  ratio_response = TRUE,
  show_ic50 = TRUE,
  add_ci = TRUE,
  transform_dose = TRUE,
  n_samples = 1000
)

## ----visualize_another, fig.height=5, fig.width=7-----------------------------
visualizeResponseProtein(
  data = dia_prepared,
  protein_name = "PROTEIN_B",
  drug_name = "Drug1",
  ratio_response = TRUE,
  show_ic50 = TRUE,
  add_ci = TRUE
)

## ----compare_scales, fig.height=4, fig.width=12-------------------------------
p1 <- visualizeResponseProtein(
  data = dia_prepared,
  protein_name = "PROTEIN_A",
  drug_name = "Drug1",
  ratio_response = FALSE,
  show_ic50 = TRUE,
  add_ci = FALSE
)

p2 <- visualizeResponseProtein(
  data = dia_prepared,
  protein_name = "PROTEIN_A",
  drug_name = "Drug1",
  ratio_response = TRUE,
  show_ic50 = TRUE,
  add_ci = FALSE
)

gridExtra::grid.arrange(p1, p2, ncol = 2)

## ----visualize_ic25, fig.height=5, fig.width=7--------------------------------
visualizeResponseProtein(
  data = dia_prepared,
  protein_name = "PROTEIN_B",
  drug_name = "Drug1",
  ratio_response = TRUE,
  show_ic50 = TRUE,
  add_ci = TRUE,
  transform_dose = TRUE,
  target_response = 0.25
)

## ----future_experiment, message=FALSE, warning=FALSE--------------------------
# pull unique doses 
dose_levels = as.numeric(unique(dia_prepared$dose))

# run simulation
custom_simulation <- futureExperimentSimulation(
  N_proteins = 300,
  N_rep = 2,
  N_Control_Rep = 2,
  Concentrations = dose_levels,
  data = dia_prepared,
  strong_proteins = "PROTEIN_B",
  weak_proteins = "PROTEIN_A",
  no_interaction_proteins = "PROTEIN_C",
  drug_name = "Drug1",
  IC50_Prediction = FALSE
)

print(custom_simulation$Hit_Rates_Plot)
print(custom_simulation$Template_Used)

## ----tpr_power_curve, eval=TRUE, message=FALSE, warning=FALSE-----------------
power_results <- run_tpr_simulation(
  rep_range = c(2,4),
  concentrations = dose_levels,
  dose_range = c(2, 9),
  data = dia_prepared,
  protein = "PROTEIN_A",
  n_proteins = 300
)

plot_tpr_power_curve(power_results)


## ----session_info-------------------------------------------------------------
sessionInfo()