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