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