Analysing a two-condition drug screen with diffHTS

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:

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()

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:

raw <- read_hts_plate(hts_primary_raw)
unique(raw$plate_id)
#> [1] "P01" "P02" "P03"

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):

layout_file <- system.file("extdata", "plate_layout_example.csv",
                           package = "diffHTS")
layout <- read_plate_layout(layout_file, plate_id = "P01")
table(layout$well_type)
#> 
#>    Blank       NC       PC compound 
#>       36        6        6       48
head(layout)
#>   plate_id well row col well_type compound_id concentration
#> 1      P01   A1   A   1     Blank        <NA>            NA
#> 2      P01   A2   A   2     Blank        <NA>            NA
#> 3      P01   A3   A   3     Blank        <NA>            NA
#> 4      P01   A4   A   4     Blank        <NA>            NA
#> 5      P01   A5   A   5     Blank        <NA>            NA
#> 6      P01   A6   A   6     Blank        <NA>            NA

# 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)
#> 
#>    Blank       NC       PC compound 
#>       36        6        6       48

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:

summarize_plate_setup(raw)
#>   plate_id n_wells plate_format n_rows n_cols n_nc n_pc n_blank n_compound
#> 1      P01      96           96      8     12    6    6      36         48
#> 2      P02      96           96      8     12    6    6      36         48
#> 3      P03      96           96      8     12    6    6      36         48
#>   n_unique_compound nc_columns pc_columns blank_columns compound_columns
#> 1                48      2, 11      2, 11          1-12             3-10
#> 2                48      2, 11      2, 11          1-12             3-10
#> 3                48      2, 11      2, 11          1-12             3-10

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

raw <- baseline_subtract(raw)
raw <- detect_outlier_wells(raw)
calc_z_prime(raw)[, c("plate_id", "z_prime", "pass")]
#>   plate_id    z_prime  pass
#> 1      P01  0.8924209  TRUE
#> 2      P02  0.8327377  TRUE
#> 3      P03 -1.0916305 FALSE
good <- filter_valid_plates(raw)
attr(good, "dropped_plates")   # P03 is dropped for failing Z'
#> [1] "P03"

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:

qc <- calc_plate_qc(raw)
qc[, c("plate_id", "cv_nc", "cv_pc", "sb", "sn",
       "z_prime", "robust_z_prime", "ssmd", "pass")]
#>   plate_id     cv_nc    cv_pc        sb       sn    z_prime robust_z_prime
#> 1      P01  2.532978 15.07000 17.717092 110.9296  0.8924209      0.9220759
#> 2      P02  4.399919 13.75264 16.442511 112.2876  0.8327377      0.8026136
#> 3      P03 28.592383 49.22758  2.892112   3.8436 -1.0916305     -0.3772639
#>        ssmd  pass
#> 1 35.313030  TRUE
#> 2 20.969913  TRUE
#> 3  1.966115 FALSE

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):

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

calc_ssmd(raw)[, c("plate_id", "ssmd", "pass")]
#>   plate_id      ssmd  pass
#> 1      P01 35.313030  TRUE
#> 2      P02 20.969913  TRUE
#> 3      P03  1.966115 FALSE
calc_cv(raw)
#>   plate_id     cv_nc    cv_pc  pass
#> 1      P01  2.532978 15.07000  TRUE
#> 2      P02  4.399919 13.75264  TRUE
#> 3      P03 28.592383 49.22758 FALSE

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:

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:

scored <- calc_well_zscore(raw)
plot_plate_qc(scored, metric = "robust_zscore")

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:

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.

norm <- norm_by_control(good)
round(range(norm$inhibition, na.rm = TRUE), 1)
#> [1]  -6.5 106.5

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.

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.

cv <- calc_replicate_cv(norm)
calc_replicate_correlation(norm)
#>   plate_x plate_y  n pearson_r r_squared
#> 1     P01     P02 48 0.9539659  0.910051

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.

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().

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.

hits <- select_primary_hit(norm, threshold = 50)
summarize_primary_hit(hits)
#>   n_compounds n_hits  hit_rate cutoff
#> 1          48     10 0.2083333     50

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:

sig <- select_sigma_hits(norm, n_sigma = 3)
head(sig)
#>   compound_id activity    sigma sigma_level       band is_hit
#> 1      CPD022 91.87347 5.073308           3 >= 3 sigma   TRUE
#> 2      CPD039 91.87205 5.073220           3 >= 3 sigma   TRUE
#> 3      CPD037 88.86450 4.885608           3 >= 3 sigma   TRUE
#> 4      CPD017 86.42806 4.733622           3 >= 3 sigma   TRUE
#> 5      CPD041 84.57373 4.617949           3 >= 3 sigma   TRUE
#> 6      CPD040 82.74610 4.503941           3 >= 3 sigma   TRUE
summarize_sigma_hits(sig)
#>   sigma n_beyond frac_beyond expected_frac expected_n
#> 1     1       16   0.3333333   0.158655254 7.61545219
#> 2     2       12   0.2500000   0.022750132 1.09200633
#> 3     3        9   0.1875000   0.001349898 0.06479511

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:

c(robust = sum(select_sigma_hits(norm, method = "robust")$is_hit),
  sd     = sum(select_sigma_hits(norm, method = "sd")$is_hit))
#> robust     sd 
#>      9      0

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:

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.

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))
#>         curve_id compound_id cell_line  n       rsq      auc      ic50
#> 1  CPD003 | A549      CPD003      A549 16 0.9837911 3.386037  6.612306
#> 2 CPD003 | H1299      CPD003     H1299 16 0.9895384 3.446830  7.703806
#> 3  CPD003 | MRC5      CPD003      MRC5 16 0.9837632 3.756423 20.884344
#> 4  CPD007 | A549      CPD007      A549 16 0.9938004 2.724255  1.126206
#> 5 CPD007 | H1299      CPD007     H1299 16 0.9893438 2.772721  1.696998
#> 6  CPD007 | MRC5      CPD007      MRC5 16 0.9913019 3.241486  3.636173
#>        emax      bottom      hill pass
#> 1 0.9917634  0.13328645 1.0436184 TRUE
#> 2 1.0059772  0.11167686 0.9192067 TRUE
#> 3 1.0217732 -0.05208097 1.0713194 TRUE
#> 4 1.0081130  0.09376054 0.9716556 TRUE
#> 5 1.0024735  0.03363457 0.9799235 TRUE
#> 6 1.0082866  0.13373237 1.3979129 TRUE

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.

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.

m  <- build_auc_matrix(drc$meta)
cl <- cluster_auc_matrix(m, k = 2)
table(cl$clusters)
#> 
#> 1 2 
#> 4 6

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.

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.

ranked <- rank_hit_compound(params)
table(ranked$tier)
#> 
#>     Weak Moderate   Strong 
#>       10       10       10
head(annotate_hit_info(ranked, hts_compound_meta))
#>   compound_id       curve_id      auc cell_line  n       rsq       ic50
#> 1      CPD037  CPD037 | A549 1.715935      A549 16 0.9961264 0.09475065
#> 2      CPD044  CPD044 | MRC5 2.685713      MRC5 16 0.9964427 0.82958180
#> 3      CPD031  CPD031 | A549 1.737770      A549 16 0.9947532 0.12298169
#> 4      CPD044  CPD044 | A549 2.044402      A549 16 0.9895371 0.28074446
#> 5      CPD037 CPD037 | H1299 2.168967     H1299 16 0.9897184 0.14304989
#> 6      CPD007  CPD007 | A549 2.724255      A549 16 0.9938004 1.12620646
#>       emax     bottom      hill pass    score   tier compound_name  target
#> 1 1.029703 0.06769499 1.0745491 TRUE 5.756930 Strong  Compound-037    mTOR
#> 2 1.043773 0.09370735 0.8004819 TRUE 4.600210 Strong  Compound-044 Unknown
#> 3 0.985824 0.06448258 1.5656243 TRUE 3.333287 Strong  Compound-031    mTOR
#> 4 1.005560 0.03771027 1.1006766 TRUE 2.589005 Strong  Compound-044 Unknown
#> 5 1.006866 0.19067325 1.1770503 TRUE 2.495076 Strong  Compound-037    mTOR
#> 6 1.008113 0.09376054 0.9716556 TRUE 2.213547 Strong  Compound-007    PARP
#>                 moa  library_source
#> 1  Kinase inhibitor Natural product
#> 2 Apoptosis inducer Natural product
#> 3    Antimetabolite  Kinase-focused
#> 4 Apoptosis inducer Natural product
#> 5  Kinase inhibitor Natural product
#> 6    Antimetabolite Natural product

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.

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.