---
title: "Analysing a two-condition drug screen with diffHTS"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Analysing a two-condition drug screen with diffHTS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(diffHTS)
```

## Overview

`diffHTS` is a one-stop toolkit for **large-scale** microplate
(96/384/1536-well) high-throughput drug screening. It is built for campaigns
that span **many plates and several experiments at once**, and covers the whole
chain from raw luminescence/fluorescence readouts through rigorous
**quality control**, replicate checks, dose-response fitting, AUC, clustering
and hit ranking. It supports both single-concentration **primary** screens and
gradient **secondary** (dose-response) screens.

Two comparisons sit at the heart of the package and shape every module:

- **Two conditions.** Drugs are compared between two experimental conditions —
  for example irradiated (`Gy2`) versus non-irradiated (`Gy0`) cells — via the
  differential (delta-AUC) workflow, to find compounds whose effect *depends* on
  the condition (e.g. radiosensitisers).
- **Cancer versus normal cell lines.** The same compounds are screened across
  cancer and normal cell lines, so that hits can be selected for a
  **therapeutic window** — potent against cancer cells while sparing normal
  cells — rather than for raw potency alone.

The analysis is organised into seven modules (each with business + companion
plotting functions), shown first as an end-to-end pipeline. The
**two-condition differential** workflow (delta-AUC between, e.g., `Gy0` vs
`Gy2`) is described afterwards.

### Publication-ready figures

Every companion plotting function returns a standard `ggplot` object drawn in a
consistent, colourblind-safe *modern academic* style: a white background with an
L-shaped axis and light dashed grid (`theme_hts()`), a low-saturation Okabe-Ito
qualitative palette (`hts_pal()`), muted semantic colours (blue for pass/NC,
vermillion for fail/PC), direct on-plot value labels and thin white bar
separators. Both helpers are exported, so you can restyle your own figures to
match:

```{r 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()
```

## The seven-module pipeline

### Module 1 — Import and pre-QC

`read_hts_plate()` ingests a plate table (data frame or csv/txt/xlsx) into a
standardised `hts_raw` object. `baseline_subtract()` removes the blank-well
background, `detect_outlier_wells()` flags anomalies by the IQR rule,
`calc_z_prime()` computes each plate's Z' factor and `filter_valid_plates()`
drops plates below the threshold.

**A screening run normally spans several plates at once, and every Module 1
function is plate-aware.** A single `hts_raw` object can hold any number of
plates stacked in its `plate_id` column; background subtraction, outlier
detection and Z' scoring are all computed *within each plate*, and
`filter_valid_plates()` drops whole failing plates while keeping the rest. The
bundled `hts_primary_raw` demonstrates this with **three 96-well plates**:

```{r m1-read}
raw <- read_hts_plate(hts_primary_raw)
unique(raw$plate_id)
```

**How the plates are set up — you define the layout.** Nothing about where the
controls, blanks and compounds sit is fixed by the package: you decide it, and
record it in a **plate-map file** you author in Excel or CSV. The map is a grid
that mirrors the physical plate — the first column holds the row letters, the
remaining headers are the plate column numbers, and each cell names what that
well contains (a negative-control label such as `DMSO`, a positive-control
label such as `kill_ctrl`, a blank/edge label such as `edgewell`, or a compound
identifier). `read_plate_layout()` turns that grid into a tidy layout, and
`apply_plate_layout()` stamps those roles onto the measured signal so the *what
each well is* and the *what was measured* stay in separate files, exactly as a
lab records them.

**Two layout rules worth following.** However you design a plate, two spatial
conventions make the data more trustworthy and are used throughout the bundled
demos:

1. **Leave the outer ring empty.** Wells on the plate edge evaporate faster than
   interior wells, which biases their readings (the "edge effect"). Reserve the
   outer ring purely as a buffer — no drugs, no controls — and fill it only with
   medium/blank. The demos leave row A, row H and columns 1 and 12 empty on the
   96-well plate (a two-well ring on 384-well).
2. **Place controls on a diagonal.** Rather than parking all negative controls
   in one column and all positive controls in another, split each control across
   two inset columns diagonally — positive controls top-left and bottom-right,
   negative controls top-right and bottom-left. This *rotationally balanced*
   arrangement averages out any left-to-right or top-to-bottom signal gradient,
   so the controls report the true plate-wide window.

The package ships a small example map that follows both rules: an empty
`edgewell` ring with diagonal controls in the two inset columns (2 and 11):

```{r 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)
```

```{r 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)
```

The bundled `hts_primary_raw` already carries its well roles, and its layout
follows the same two rules: the outer ring (row A, row H, columns 1 and 12) is
an empty `Blank` evaporation buffer, the two inset columns 2 and 11 hold the
controls on a diagonal, and the 48 compounds fill the interior (rows B–G,
columns 3–10). `P01` and `P02` are good replicate plates; `P03` was
deliberately simulated with poor control separation so that QC can reject it.
`summarize_plate_setup()` reads whatever layout each plate actually uses
straight off the data — one row per plate — so you can confirm every plate was
built the way you intended before going further:

```{r m1-setup}
summarize_plate_setup(raw)
```

With the layout confirmed, run the per-plate pre-QC over all three plates at
once:

```{r 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'
```

**A comprehensive quality panel.** The Z'-factor is only one of several
complementary ways to judge a plate. `calc_plate_qc()` returns the whole panel
in one per-plate table — the classic **Z'-factor** and its outlier-resistant
**robust Z'-factor** (median/MAD based), the sample-aware **Z-factor**,
**SSMD**, the **signal-to-background (S/B)** and **signal-to-noise (S/N)**
ratios, and the **CV%** of each control:

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

Each metric answers a different question, so a plate can pass one and fail
another (here `P03` keeps an adequate S/B dynamic range yet fails everything
that accounts for variability):

| Metric | Function | What it tells you | Good |
|---|---|---|---|
| Z'-factor | `calc_z_prime()` | Overall control-only assay quality | >= 0.5 |
| Robust Z'-factor | `calc_robust_z_prime()` | Same, but immune to outlier control wells | >= 0.5 |
| Z-factor | `calc_z_factor()` | Window vs. the actual library spread (often low in primary screens) | >= 0.5 |
| SSMD | `calc_ssmd()` | Effect-size of control separation, for hit cut-offs | abs >= 2 |
| S/B | `calc_sb_ratio()` | Raw dynamic range of the assay | >= 2 |
| S/N | `calc_sn_ratio()` | Assay window relative to background noise | >= 3 |
| CV% | `calc_cv()` | Control reproducibility (pipetting/reagents) | <= 20% |

Each metric also has its own standalone function (mirroring `calc_z_prime()`)
for when you need just one number with its own threshold:

```{r m1-qc-single}
calc_ssmd(raw)[, c("plate_id", "ssmd", "pass")]
calc_cv(raw)
```

`plot_plate_qc()` is a single visualiser: pick any metric with the `metric`
argument and it draws the plates with the right acceptance line and pass/fail
colouring. The failing plate `P03` stands out on every metric:

```{r 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
```

Per-well **Z-scores** and **robust Z-scores** (the median/MAD score preferred in
HTS because hits do not distort it) come from `calc_well_zscore()` and plug into
the same visualiser:

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

```{r m1-heatmap, fig.width = 8, fig.height = 5}
plot_plate_heatmap_raw(raw, plate = "P01")
```

The negative (`NC`) and positive (`PC`) control wells are automatically circled
so they can be told apart from the drug-treated wells. The same works for
larger formats. The bundled `hts_plate_384` dataset mirrors a real 384-well
screen (EXP87) and applies exactly the same two rules: the outer two-well ring
(rows A, B, O, P and columns 1, 2, 23, 24) is kept free of compound to avoid
edge-evaporation artefacts. Those wells still hold medium but no live cells, so
they read at the same low, "dead" level as the kill control (dark, not missing).
The controls sit on a diagonal in the two inset
columns 3 and 22 — column 3 holds the positive (kill) control in its top half
(rows C–H) and the negative (vehicle) control in its bottom half (rows I–N),
while column 22 mirrors it. The circled controls trace out the tell-tale
diagonal:

```{r 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)
```

### Module 2 — Normalisation and activity

Raw signals cannot be compared across plates: each plate has its own cell
number, reagent batch and instrument gain. `norm_by_control()` removes that
plate-to-plate offset by rescaling every well against **that plate's own
controls**, so all plates land on a common 0–100% activity scale. With the
negative control (vehicle, full growth) defining 0% inhibition and the positive
control (kill) defining 100% inhibition, a well with signal `x` is scored as

`inhibition (%) = 100 * (mean(NC) - x) / (mean(NC) - mean(PC))`

so a well behaving like the vehicle scores about 0, a fully killed well about
100, and viability is simply `100 - inhibition`. Because the controls are taken
per plate, this both normalises the scale **and** corrects for plate-level
effects in one step. `merge_plate_data()` then stacks the normalised plates into
a single library-wide table.

```{r m2}
norm <- norm_by_control(good)
round(range(norm$inhibition, na.rm = TRUE), 1)
```

**Why plot the inhibition histogram?** `plot_inhibition_hist()` shows the
distribution of per-well inhibition across a plate, and it is the quickest sanity
check that normalisation worked. In a well-behaved primary screen the vast
majority of compounds are inactive, so the histogram should show a **tall peak
near 0% inhibition** (the inactive library) with only a **thin right-hand tail**
of active compounds approaching 100%. A peak sitting away from 0, a bimodal
shape, or a shifted centre signals a normalisation or control problem on that
plate before any hit is called. `plot_plate_heatmap_inhibition()` shows the same
values back in their physical well positions to expose any spatial pattern.

```{r m2-plots, fig.width = 7, fig.height = 4}
plot_plate_heatmap_inhibition(norm, plate = "P01")
plot_inhibition_hist(norm, plate = "P01")
```

### Module 3 — Replicate consistency

Screening decisions are only as trustworthy as the measurements behind them, so
before calling any hits we check that independent replicates of the same
compound agree. `calc_replicate_cv()` computes each compound's **coefficient of
variation** (CV = SD / mean) across its replicate wells, and
`calc_replicate_correlation()` reports the plate-to-plate Pearson correlation;
`filter_bad_replicate()` then removes compounds whose replicates disagree.

```{r m3}
cv <- calc_replicate_cv(norm)
calc_replicate_correlation(norm)
```

`plot_replicate_scatter()` draws every pairwise replicate comparison at once
(P01 vs P02, P01 vs P03, P02 vs P03) in one faceted figure. Following the house
style of the screening reports, points are coloured by well type (compounds
blue, negative controls black, positive controls red, blanks gold), each panel
carries a regression line with a shaded 95% confidence band, and the Pearson
`R` and its p-value are annotated. Pass `plate_x` and `plate_y` to focus on a
single pair.

```{r m3-plots, fig.width = 9, fig.height = 3.4}
plot_replicate_scatter(norm)
```

**Why look at the CV distribution?** A single global CV cut-off hides how the
whole library behaves. Plotting the distribution of per-compound CVs shows
whether the assay is tight (most mass well below the threshold) or noisy (a long
tail crossing it), reveals whether a small number of erratic compounds — rather
than a systemic problem — drive the failures, and lets you set an
evidence-based, rather than arbitrary, acceptance threshold (dashed line, `15%`
by convention). Compounds to the right of the line are dropped by
`filter_bad_replicate()`.

```{r m3-cv, fig.width = 6, fig.height = 4}
plot_cv_distribution(cv, cv_max = 15)
```

### Module 4 — Primary hit selection

`select_primary_hit()` flags actives by a fixed threshold or a percentile rank;
`summarize_primary_hit()` reports the hit rate.

```{r m4}
hits <- select_primary_hit(norm, threshold = 50)
summarize_primary_hit(hits)
```

```{r m4-plots, fig.width = 6, fig.height = 4}
plot_primary_inhibition_rank(hits)
plot_hit_bar_count(hits)
```

`plot_hit_bar_count()` simply tallies **hit vs non-hit** from any table carrying
an `is_hit` column, so it works for *both* ways of calling hits: the
**rank/threshold** method above (`select_primary_hit()`) and the
**distribution (sigma)** method below (`select_sigma_hits()`). Drawing it for
each makes the two definitions directly comparable — a fixed threshold and a
sigma tail usually agree on the strongest actives but disagree at the margin, so
seeing both counts side by side tells you how sensitive the hit list is to the
selection rule.

**A distribution-based (sigma) alternative.** A fixed cut-off ignores how noisy
the library actually is. The classic HTS view instead treats the per-compound
activity of the whole library as an approximately normal *inactive*
distribution and calls anything far out in the tail — conventionally beyond
**3 sigma** — a hit. `select_sigma_hits()` scores every compound by how many
standard deviations it sits from the library centre and bins it by whether it
clears 1, 2 or 3 sigma:

```{r m4-sigma}
sig <- select_sigma_hits(norm, n_sigma = 3)
head(sig)
summarize_sigma_hits(sig)
```

The `summarize_sigma_hits()` table contrasts the observed tail counts with what
a perfectly normal (hit-free) library would predict, so the enrichment beyond
3 sigma is obvious — far more compounds sit in the tail than chance allows.

By default the centre and spread are estimated **robustly** (median and MAD),
which is deliberate. A handful of genuine hits inflate a plain mean/SD so much
that they hide themselves: with `method = "sd"` the standard deviation is
stretched by the very hits we are hunting and *nothing* reaches 3 sigma, whereas
the outlier-resistant median/MAD keeps a tight null and recovers the real tail:

```{r m4-sigma-compare}
c(robust = sum(select_sigma_hits(norm, method = "robust")$is_hit),
  sd     = sum(select_sigma_hits(norm, method = "sd")$is_hit))
```

`plot_sigma_hits()` shows the whole picture: the activity histogram coloured by
sigma band, the fitted normal curve, and the 1/2/3 sigma cut-offs. The bump of
orange/red bars poking out past the normal curve on the right is the hit tail.
Passing the sigma table to `plot_hit_bar_count()` gives the matching hit vs
non-hit tally for this distribution-based rule, directly comparable to the
threshold-based count above:

```{r m4-sigma-plot, fig.width = 7, fig.height = 4.5}
plot_sigma_hits(sig)
plot_hit_bar_count(sig)
```

### Module 5 — Dose-response fitting and AUC

A primary screen at a single concentration only tells you *whether* a compound is
active; to know *how* active, the confirmed hits are re-tested across a
concentration gradient and fitted with a dose-response curve. `import_dose_response()`
standardises the gradient data, `fit_4pl_curve()` fits a 4PL model per compound
(and group, here cell line), `extract_drc_params()` pulls IC50/Emax/Hill and
`calc_drc_auc()` integrates each curve.

**Why compute both IC50 and AUC?** They describe *different* aspects of a
compound and neither is sufficient alone:

- **IC50** (the concentration giving half-maximal effect) measures **potency** —
  how little compound is needed. It is read from a single point on the curve and
  ignores what happens at the plateau, so a very potent compound that only kills
  40% of cells still gets a small IC50.
- **Emax / Hill slope** describe the **efficacy** (maximum achievable effect) and
  the steepness of the response.
- **AUC** (area under the curve) integrates the whole curve, combining potency
  *and* efficacy into one robust number. Because it uses every concentration it
  is far less sensitive to a single noisy point or to a poorly-constrained IC50
  when a curve does not fully plateau, which makes it the more stable ranking
  metric in HTS.

```{r 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))
```

```{r 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])
```

Because the 4PL is fitted **per cell line**, `plot_batch_drc_overlay()` lets you
overlay the *same* compound's curves across the different cell lines on one axis.
This is the quickest way to read selectivity directly off the dose-response: a
curve that drops steeply in the cancer lines but stays high (flat, near full
viability) in the normal line marks a compound that kills cancer cells while
sparing normal ones — exactly the profile the screen is looking for.

**Why correlate IC50 against AUC?** Because the two metrics capture potency and
integrated activity, comparing them is a built-in consistency check. For most
compounds a lower IC50 (more potent) goes with a larger integrated effect, so
the two are correlated — `plot_ic50_auc_cor()` draws the trend line and prints
the Pearson `R` and p-value (computed on `log10(IC50)` to match the log x-axis)
directly on the figure. The value of the plot is in the **outliers**: a compound
that is potent (small IC50) but sits off the trend with a modest AUC is efficacy-
limited (steep but shallow curve), while a weak-IC50 compound with a
surprisingly large AUC has a broad, gradual response. These are exactly the
molecules worth a second look, and the plot separates strong / moderate / weak
candidates at a glance.

```{r 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)
```

### Module 6 — AUC matrix and clustering

`build_auc_matrix()` reshapes AUCs into a compound-by-sample matrix,
`cluster_auc_matrix()` clusters it and `extract_cluster_hit()` pulls the most
active cluster.

```{r m6}
m  <- build_auc_matrix(drc$meta)
cl <- cluster_auc_matrix(m, k = 2)
table(cl$clusters)
```

Reshaping every compound's AUC into one **compound-by-cell-line matrix** is what
turns selectivity into a *batch* read-out. Each row is a compound, each column a
cell line, and clustering the rows groups compounds by their **whole
cross-cell-line profile** at once — you no longer inspect drugs one at a time.
In the row-scaled heatmap below, a compound that is active only in the cancer
columns (dark) while staying inactive in the normal `MRC5` column (light) stands
out immediately as a differential, cancer-selective hit; broadly toxic compounds
are dark across *all* columns and are easy to set aside. That side-by-side
contrast across the full library is what lets a large screen be triaged quickly.

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

### Module 7 — Ranking, annotation and report

`rank_hit_compound()` combines the metrics into a Strong/Moderate/Weak score,
`annotate_hit_info()` joins compound metadata, `export_hit_table()` writes a
CSV, and `generate_hts_report()` renders a full HTML report.

```{r m7}
ranked <- rank_hit_compound(params)
table(ranked$tier)
head(annotate_hit_info(ranked, hts_compound_meta))
```

A trustworthy final hit is not defined by any single test but by **agreement
across the evidence gathered in the previous modules**:

1. **Primary rank / threshold** (Module 4) — the compound is active in the
   single-concentration screen.
2. **Sigma tail** (Module 4) — its activity sits well beyond the library's null
   distribution (e.g. past 3 sigma), so it is unlikely to be noise.
3. **Cancer-versus-normal differential** (Modules 5–6 and the two-condition
   workflow below) — it is potent in the cancer lines/treated condition while
   sparing the normal line, i.e. it has a therapeutic window.

`rank_hit_compound()` produces the tiered score, `annotate_hit_info()` attaches
the library metadata, and the compounds that clear **all three** criteria are
the ones worth carrying forward — a shortlist that a large primary library can
be reduced to with confidence.

## Two-condition differential workflow

The remainder of this vignette covers the original differential (delta-AUC)
analysis used to compare drug sensitivity between **two conditions** -- for
example irradiated (`Gy2`) versus non-irradiated (`Gy0`) cells, or a cancer
versus a normal cell line.

### Experimental design: cancer versus normal cell lines

The screen is built around a **therapeutic-window** question rather than raw
potency. Each compound is tested in parallel on one or more **cancer** cell
lines (here `A549`, `H1299`, `H1975`, `H2030`) and on a **normal** cell line
(here `MRC5`), and the same 4PL dose-response fit gives an AUC (or IC50) per
compound per line. A useful hit is one that is **potent against the cancer
cells** (low cancer AUC/IC50, i.e. it kills them) while **sparing the normal
cells** (high `MRC5` AUC/IC50). Comparing the two rather than looking at a
single line filters out compounds that are simply broadly cytotoxic: those lower
the AUC everywhere, including in `MRC5`, and offer no selectivity. In the demo
data `MRC5` is deliberately made more resistant (a higher IC50) than the cancer
lines, so selective compounds separate out as a large gap between the cancer and
normal AUC/IC50. The same logic applies to the radiation setting: keep drugs
that sensitise the treated condition (`Gy2`) while leaving the untreated
baseline (`Gy0`) largely unaffected.

## 1. Dose-response and AUC

`four_pl()` evaluates the four-parameter logistic model and `compute_auc()`
integrates it. `fit_dose_response()` fits a curve directly from well-level data
using a self-contained base-R optimiser (no external fitting package required),
returning a `dsr_4pl` object with the usual `coef()`, `predict()`, `fitted()`
and `residuals()` methods.

```{r 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)
```

## 2. Quality control

`calculate_qc_metrics()` derives the standard plate metrics (Z-factor, Z',
signal-to-background, signal-to-noise, SSMD) from per-control summary
statistics, and `flag_qc_plates()` applies acceptance thresholds.

```{r 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)
```

```{r qc-boxplot}
plot_qc_boxplot(screen_doseresponse)
```

`plot_plate_heatmap()` shows a whole plate in its physical layout, colouring
each well by cell growth. The demo plate follows standard HTS practice, so three
zones are visible at a glance: an **empty evaporation-buffer ring** around the
outer edge (row A, row H, columns 1 and 12), which holds no live cells and so
reads dark like the kill control; the **controls** placed in the two flanking
columns (2 and 11) in a rotationally balanced *diagonal* pattern -- positive
(kill) controls top-left and bottom-right, negative (DMSO) controls top-right
and bottom-left -- so any spatial gradient averages out; and the **dose series**
of one drug per row filling the interior block (rows B-G, columns 3-10). Keeping
compounds off the edge avoids the well-known edge effect (faster evaporation and
uneven readouts around the plate rim).

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

For a whole screen of many plates, `screen_plate_qc` holds one row per plate
with its control readouts and QC metrics. `plot_control_scatter()` compares the
negative- and positive-control readouts across plates, and
`plot_qc_metric_trend()` tracks any single metric (here Z') with failing plates
highlighted.

```{r 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)
```

When several experiments are screened together, circular summaries make the
whole campaign legible at once. `plot_qc_metric_circular()` places one sector
per experiment and one thin bar per plate rising from a zero baseline to its
Z' value, with the value printed at each bar tip (vermillion = fails the Z'
threshold, blue = passes) and a generous outer margin so no bar touches the
sector edge, while `plot_qc_circular_heatmap()` shows every control well of
every plate as a ring.

```{r 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)
```

## 3. Differential AUC between the two conditions

Given a long table of per-drug, per-condition AUCs, `compute_delta_auc()`
returns one row per drug with `delta_auc` (`treatment - baseline`) and its
z-score.

```{r 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"))
```

## 4. Hit selection

Across several cell lines, `select_hits_cutoff()` keeps drugs below a threshold
in enough cell lines, while `select_hits_sigma()` keeps drugs beyond a sigma
cut-off. Here the cancer lines are scored together: a hit must have a negative
delta-AUC (more sensitive under treatment) in at least four of them. Reading the
result alongside the normal line `MRC5` is what reveals selectivity -- prefer
compounds that pass in the cancer lines but leave `MRC5` near zero (spared).

```{r 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)
```

## 5. Visualisation

`plot_dose_response_curves()` overlays the observed responses and fitted 4PL
curves for a single drug under both conditions. In the demo data `DrugA` is a
strong radiosensitiser — its `Gy2` IC50 is about ten-fold lower than its `Gy0`
IC50 — so the two curves are **clearly separated**: the treated (`Gy2`) curve
sits well to the left, killing cells at much lower concentrations than the
untreated (`Gy0`) curve. That visible left-shift is the signature of a compound
whose effect depends on the condition:

```{r curves}
drugA <- subset(
  screen_doseresponse,
  drug_name == "DrugA" & experiment_type == "sample"
)
plot_dose_response_curves(drugA, drug = "DrugA")
```

`plot_condition_scatter()` compares each drug across two cell lines at once, with
the axes centred on the origin so the **four quadrants are equal in size** and
easy to read. Here the normal line `MRC5` is on the x-axis and the cancer line
`A549` on the y-axis, and the values are delta-AUCs (negative = more killing).
The quadrant a drug falls in tells you its selectivity at a glance:

- **Lower-right** (normal `MRC5` change near zero or positive = *spared*, cancer
  `A549` change strongly negative = *killed*): the therapeutic-window drugs — the
  ones the screen is after.
- **Lower-left**: kills both cancer and normal cells (broadly toxic,
  non-selective).
- **Upper-right / near the origin**: little effect on either line.

```{r scatter, fig.width = 5.5, fig.height = 5.5}
plot_condition_scatter(screen_delta_auc, x = "MRC5", y = "A549")
```

Finally, `plot_delta_auc_heatmap()` lays the whole library out as a
compound-by-cell-line grid with the cancer and normal columns annotated
separately. Cancer-selective hits appear as rows that are **strongly coloured
across the cancer columns but pale in the `MRC5` (normal) column**, so the
differential between normal and cancer cells is visible for every compound at
once — turning a large screen into a quick, side-by-side triage (requires the
suggested `ComplexHeatmap` and `circlize` packages):

```{r 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")
)
```
