## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6.5,
  fig.height = 4,
  dpi = 120,
  dev = "png"
)

## ----setup--------------------------------------------------------------------
library(diffHTS)

## ----styling------------------------------------------------------------------
library(ggplot2)
ggplot(mtcars, aes(factor(cyl), mpg, fill = factor(cyl))) +
  geom_boxplot(color = "white", linewidth = 0.3) +
  scale_fill_manual(values = hts_pal()) +
  labs(x = "Cylinders", y = "MPG", fill = "cyl") +
  theme_hts()

## ----m1-read------------------------------------------------------------------
raw <- read_hts_plate(hts_primary_raw)
unique(raw$plate_id)

## ----m1-layout----------------------------------------------------------------
layout_file <- system.file("extdata", "plate_layout_example.csv",
                           package = "diffHTS")
layout <- read_plate_layout(layout_file, plate_id = "P01")
table(layout$well_type)
head(layout)

## ----m1-apply-----------------------------------------------------------------
# A signal export carries only well + reading; the layout supplies the roles.
set.seed(1)
signal_export <- data.frame(well = layout$well,
                            signal = runif(nrow(layout), 1e4, 1e5))
raw_mapped <- apply_plate_layout(signal_export, layout, plate_id = "P01")
table(raw_mapped$well_type)

## ----m1-setup-----------------------------------------------------------------
summarize_plate_setup(raw)

## ----m1-----------------------------------------------------------------------
raw <- baseline_subtract(raw)
raw <- detect_outlier_wells(raw)
calc_z_prime(raw)[, c("plate_id", "z_prime", "pass")]
good <- filter_valid_plates(raw)
attr(good, "dropped_plates")   # P03 is dropped for failing Z'

## ----m1-qc--------------------------------------------------------------------
qc <- calc_plate_qc(raw)
qc[, c("plate_id", "cv_nc", "cv_pc", "sb", "sn",
       "z_prime", "robust_z_prime", "ssmd", "pass")]

## ----m1-qc-single-------------------------------------------------------------
calc_ssmd(raw)[, c("plate_id", "ssmd", "pass")]
calc_cv(raw)

## ----m1-qc-plots, fig.width = 6, fig.height = 3.5-----------------------------
plot_plate_qc(qc, metric = "z_prime")        # assay window (controls only)
plot_plate_qc(qc, metric = "robust_z_prime") # robust to outlier control wells
plot_plate_qc(qc, metric = "ssmd")           # control separation
plot_plate_qc(qc, metric = "sn")             # signal-to-noise
plot_plate_qc(qc, metric = "cv")             # control CV%, both controls

## ----m1-zscore, fig.width = 6, fig.height = 3.5-------------------------------
scored <- calc_well_zscore(raw)
plot_plate_qc(scored, metric = "robust_zscore")

## ----m1-heatmap, fig.width = 8, fig.height = 5--------------------------------
plot_plate_heatmap_raw(raw, plate = "P01")

## ----m1-384, fig.width = 9, fig.height = 5------------------------------------
plot_plate_heatmap(hts_plate_384, fill = "signal", well_col = "well",
  control_col = "well_type", plate_col = "plate_id", plate_type = 384)

## ----m2-----------------------------------------------------------------------
norm <- norm_by_control(good)
round(range(norm$inhibition, na.rm = TRUE), 1)

## ----m2-plots, fig.width = 7, fig.height = 4----------------------------------
plot_plate_heatmap_inhibition(norm, plate = "P01")
plot_inhibition_hist(norm, plate = "P01")

## ----m3-----------------------------------------------------------------------
cv <- calc_replicate_cv(norm)
calc_replicate_correlation(norm)

## ----m3-plots, fig.width = 9, fig.height = 3.4--------------------------------
plot_replicate_scatter(norm)

## ----m3-cv, fig.width = 6, fig.height = 4-------------------------------------
plot_cv_distribution(cv, cv_max = 15)

## ----m4-----------------------------------------------------------------------
hits <- select_primary_hit(norm, threshold = 50)
summarize_primary_hit(hits)

## ----m4-plots, fig.width = 6, fig.height = 4----------------------------------
plot_primary_inhibition_rank(hits)
plot_hit_bar_count(hits)

## ----m4-sigma-----------------------------------------------------------------
sig <- select_sigma_hits(norm, n_sigma = 3)
head(sig)
summarize_sigma_hits(sig)

## ----m4-sigma-compare---------------------------------------------------------
c(robust = sum(select_sigma_hits(norm, method = "robust")$is_hit),
  sd     = sum(select_sigma_hits(norm, method = "sd")$is_hit))

## ----m4-sigma-plot, fig.width = 7, fig.height = 4.5---------------------------
plot_sigma_hits(sig)
plot_hit_bar_count(sig)

## ----m5-----------------------------------------------------------------------
dr  <- import_dose_response(hts_dose_response, response_col = "viability",
                            group_cols = "cell_line")
drc <- fit_4pl_curve(dr, group_cols = "cell_line")
drc <- calc_drc_auc(drc)
head(extract_drc_params(drc))

## ----m5-drc, fig.width = 6, fig.height = 4------------------------------------
plot_single_drc(drc, drc$meta$curve_id[1])
plot_batch_drc_overlay(drc, drc$meta$curve_id[1:4])

## ----m5-cor, fig.width = 6, fig.height = 4------------------------------------
params <- merge(extract_drc_params(drc), drc$meta[, c("curve_id", "auc")])
plot_ic50_auc_cor(params)

## ----m6-----------------------------------------------------------------------
m  <- build_auc_matrix(drc$meta)
cl <- cluster_auc_matrix(m, k = 2)
table(cl$clusters)

## ----m6-plots, fig.width = 6, fig.height = 5, eval = requireNamespace("ComplexHeatmap", quietly = TRUE)----
plot_auc_heatmap(build_auc_matrix(drc$meta, zscore = "row"))

## ----m7-----------------------------------------------------------------------
ranked <- rank_hit_compound(params)
table(ranked$tier)
head(annotate_hit_info(ranked, hts_compound_meta))

## ----doseresponse-------------------------------------------------------------
four_pl(x = c(0.1, 1, 10), top = 1, ic50 = 1, hill = 1, bottom = 0)

compute_auc(
  min_concentration = 0.01, max_concentration = 100,
  top = 1, ic50 = 5, hill = 1, bottom = 0.05
)

one <- subset(
  screen_doseresponse,
  drug_name == "DrugA" & condition == "Gy0" & experiment_type == "sample"
)
fit <- fit_dose_response(one, "normalized_cell_count", "concentration")
coef(fit)

## ----qc-----------------------------------------------------------------------
head(screen_doseresponse)

plate <- data.frame(
  experiment_type = c("negative_ctrl", "positive_ctrl", "sample"),
  mean = c(1.00, 0.05, 0.55),
  sd = c(0.05, 0.02, 0.20)
)
calculate_qc_metrics(plate)

## ----qc-boxplot---------------------------------------------------------------
plot_qc_boxplot(screen_doseresponse)

## ----plate-map, fig.width = 8, fig.height = 5---------------------------------
plot_plate_heatmap(subset(screen_plate_layout, plateID == "EXP99_Gy0"),
  control_col = "experiment_type")

## ----qc-controls, fig.width = 7, fig.height = 4-------------------------------
plot_control_scatter(screen_plate_qc)

plot_qc_metric_trend(screen_plate_qc, metric = "z_prime", threshold = 0.5)

## ----qc-circular, fig.width = 6, fig.height = 6, eval = requireNamespace("circlize", quietly = TRUE)----
plot_qc_metric_circular(screen_plate_qc, metric = "z_prime")

plot_qc_circular_heatmap(screen_plate_qc)

## ----delta--------------------------------------------------------------------
auc_long <- data.frame(
  drug_name = rep(c("DrugA", "DrugB", "DrugC", "DrugD"), each = 2),
  condition = rep(c("Gy0", "Gy2"), times = 4),
  auc = c(5.0, 2.0, 4.2, 4.1, 6.0, 3.4, 3.0, 3.1)
)
compute_delta_auc(auc_long, conditions = c("Gy0", "Gy2"))

## ----hits---------------------------------------------------------------------
head(screen_delta_auc)

cancer <- c("A549", "H1299", "H1975", "H2030")
select_hits_cutoff(screen_delta_auc, score_cols = cancer,
                   cutoff = 0, min_pass = 4)

## ----curves-------------------------------------------------------------------
drugA <- subset(
  screen_doseresponse,
  drug_name == "DrugA" & experiment_type == "sample"
)
plot_dose_response_curves(drugA, drug = "DrugA")

## ----scatter, fig.width = 5.5, fig.height = 5.5-------------------------------
plot_condition_scatter(screen_delta_auc, x = "MRC5", y = "A549")

## ----heatmap, eval = requireNamespace("ComplexHeatmap", quietly = TRUE) && requireNamespace("circlize", quietly = TRUE)----
plot_delta_auc_heatmap(
  screen_delta_auc,
  value_cols = c("A549", "H1299", "H1975", "H2030", "MRC5"),
  cell_types = c(A549 = "Cancer", H1299 = "Cancer", H1975 = "Cancer",
                 H2030 = "Cancer", MRC5 = "Normal")
)