---
title: "End to End: From Segmented Exports to a Finished Report"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{End to End: From Segmented Exports to a Finished Report}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
resource_files:
- figures/logo.svg
- figures/logo-dark.svg
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE, comment = "#>",
fig.width = 7, fig.height = 4.2, fig.align = "center",
dpi = 96
)
options(width = 84)
library(cellreportR)
```
`vignette("getting-started")` runs the short path on a single plate. This
vignette runs the long one, in the order the package is meant to be used:
```
files on disk -> units -> quality control -> standardization
-> effect sizes at two levels -> sample size -> report
```
Every stage is a separate exported function. Nothing here is a wrapper that
hides a decision, because each of these stages *is* a decision, and the
package's job is to make it visible and auditable afterwards.
All data below is synthetic and generated in-session under an explicit seed.
## 1. Ingest: the design lives in the tree and the file name
Segmentation software writes one export per acquisition, and the design facts
that matter — which compound, which experiment, which plate, which
pre-treatment interval, which exposure level — are usually encoded in the
directory layout and the file name rather than inside the file.
`cr_example_exports()` writes such a tree so the ingest functions can be
demonstrated on real files.
```{r write-tree}
root <- file.path(tempdir(), "cr_vignette_exports")
files <- cr_example_exports(root, seed = 42, n_cells = 60)
basename(files)
```
Three things in those names are not decoration:
* `(split)` marks the second pass of a **two-pass acquisition**. Left
unmerged, that unit contributes twice the cells of its neighbours.
* `(repeat)` marks a **repeated read** of a unit that already has a plain
sibling file.
* `_2.2` *looks* like a second pass of unit `2`, but is a **different physical
unit** on a different plate. It must not merge.
* `(no reagent)` marks the **specificity arm**, where the detection reagent
was left out. Those acquisitions are never samples.
### Column names, path levels, file-name grammar
Three small declarations describe the tree. Each is an object you can print
and check before a single file is read.
```{r spec}
map <- cr_column_map(
exact = c("Event Label" = "cell_id",
"Position X [px]" = "x",
"Position Y [px]" = "y",
"Nuclei - Signal Mean" = "nuclear_signal",
"Target - Signal Mean" = "target_signal"),
prefix = c("^Nuclei - Area" = "area")
)
grammar <- cr_filename_grammar(
tokens = list(interval = "[0-9]+min",
dose = "[0-9]+uM",
mode = "treated|vehicle"),
defaults = list(dose = "0uM", mode = "vehicle"),
core_patterns = c("^[0-9]+min_vehicle$", "^[0-9]+min_[0-9]+uM_treated$"),
prefix_strip = c("CompoundA", "CompoundB")
)
markers <- cr_marker_rules(
merge_unit = "\\(split\\)",
partial_plate = "\\(partial\\)",
omitted_reagent = "\\(no reagent\\)",
reacquisition = "\\(repeat\\)"
)
spec <- cr_path_spec(
levels = c("run", "compound", "experiment", "plate"),
grammar = grammar,
markers = markers
)
spec
```
`core_patterns` is the part worth pausing on. The *absence* of a token is
meaningful here: a name with no exposure token is a vehicle control. A typo in
a token therefore matches nothing, falls through to the default, and silently
reclassifies a treated acquisition as a control — pulling it into the control
denominator of its own batch. The whitelist turns that silent
misclassification into a hard error.
```{r read}
cells <- cr_read_exports(root, column_map = map, spec = spec, progress = FALSE)
dim(cells)
all(cells$parse_ok)
```
Each file now carries the design facts recovered from its path and its
name. The tokens:
```{r tokens}
per_file <- unique(cells[, c("source_file", "compound", "plate", "interval",
"dose", "mode", "replicate", "merge_unit",
"reacquisition", "omitted_reagent")])
per_file$source_file <- sub("\\.csv$", "", per_file$source_file)
per_file[, c("source_file", "interval", "dose", "mode", "replicate")]
```
and the markers, each one a typed flag rather than a string somebody has to
remember:
```{r markers-table}
per_file[, c("source_file", "merge_unit", "reacquisition", "omitted_reagent")]
```
### Files are not units
The specificity arm is split off first, so it cannot be merged into a sample
unit by accident. What remains is resolved into analysis units.
```{r units}
aside <- cells[cells$omitted_reagent, , drop = FALSE]
samples <- cells[!cells$omitted_reagent, , drop = FALSE]
samples <- cr_assign_units(
samples,
key_vars = c("compound", "experiment", "plate", "interval", "dose", "mode"),
rules = cr_merge_rules()
)
attr(samples, "n_units")
```
`cr_unit_map()` is the receipt: which unit came from how many files.
```{r unit-map}
umap <- cr_unit_map(samples)
umap$file <- sub("\\.csv$", "", basename(umap$source_path))
umap[, c("file", "n_cells", "n_files", "merged")]
length(unique(umap$well_id))
```
Nine sample files became seven units. The `(split)` and `(repeat)` files
merged into their siblings; `_2.2` stayed a unit of its own.
```{r dataset}
ds_files <- cr_dataset(samples, unit_var = "well_id")
ds_files
```
## 2. Scaling up
That tree has one unit per arm, which is enough to demonstrate ingest and not
enough to estimate anything. The rest of this vignette uses
`cr_example_screen()`, which produces the same object shape — units carrying
file provenance, a vehicle control inside every batch, a set-aside specificity
arm — at a size that supports estimation.
```{r screen}
screen <- cr_example_screen(seed = 1, n_compounds = 4, n_cells_per_well = 25)
screen
batch <- screen$batch_vars
batch
```
A batch here is the *combination* of four columns, not one column. Everything
downstream that references a control references the control of the unit's own
batch.
```{r design-table}
table(screen$design$compound, screen$design$treatment)
```
## 3. Quality control at the cell level
Sub-threshold debris first, then balancing — in that order, because excluding
first changes the pool that balancing samples from.
```{r qc-cells}
ds <- cr_exclude_small(screen, var = "area", probs = 0.10)
ds <- cr_balance_cells(ds, n_max = 30, seed = 1)
cr_qc_summary(ds)[, c("step", "cells_before", "cells_after", "percent_removed")]
```
The QC log is append-only and travels with the object, so the disposition of
every cell is recoverable at the end from the object itself rather than from
a lab notebook.
## 4. Standardization against each unit's own control
The second experiment ran at a higher raw baseline than the first, and the
second plate higher than the first. Both offsets act on control and treated
units alike, which is exactly what standardizing against the control of a
unit's own batch removes.
```{r batch-ref}
ref <- cr_batch_reference(ds, channel = "target_signal",
control_level = "Vehicle", batch_vars = batch)
head(ref[, c("batch_key", "n_cells", "ctrl_n", "ctrl_median", "has_control")], 4)
```
Inspect this table before standardizing, not after: a batch with very few
control cells, or none at all, is visible here while it is still cheap to
notice.
```{r standardize}
ds <- cr_standardize_batch(ds,
channel = "target_signal",
control_level = "Vehicle",
batch_vars = batch,
method = "log2_fc",
value_to = "log2_fc")
round(summary(ds$cells$log2_fc), 3)
```
`log2_fc` is now a per-cell value on a common scale: zero is the control of
that cell's own batch.
## 5. The unit-level gate, and what it costs
A control-referenced gate is a biological rather than a statistical
criterion: a treated unit should carry more target signal than the vehicle
control acquired in its own batch.
```{r gate}
gate <- cr_qc_gate(ds,
channel = "target_signal",
control_level = "Vehicle",
batch_vars = batch)
gate
```
The gate can only ever drop *low* units. A gate that is too strict therefore
manufactures apparent treatment effects rather than removing artefacts, so
the leverage of the exclusions has to be quantified **before** they are
applied, while both states are still in hand.
```{r gate-impact}
impact <- cr_qc_gate_impact(ds, gate,
value = "log2_fc",
group_var = "treatment",
reference_level = "Vehicle",
by = "compound")
# a narrow view of the same table: every contrast is against "Vehicle",
# so that part of the label carries no information here
data.frame(
compound = impact$compound,
contrast = sub("^Vehicle -> ", "", impact$contrast),
n_excluded = impact$n_units_excluded,
d_with = round(impact$estimate_with_excluded, 2),
d_without = round(impact$estimate_without_excluded, 2),
sign_changed = impact$sign_changed
)
```
`d_with` is the effect estimated with the failing units retained, `d_without`
the same estimate once they are dropped. Read the last column. For the three responsive compounds the gate removes a
few units and moves the estimate in the direction it already pointed. For the
weakest compound the gate removes most of the treated units and *flips the
sign*: there, the gate is doing the selecting, not the compound.
That is a reviewed decision, not an automatic one, and `cr_apply_gate()`
accepts a reviewed exclusion list rather than insisting on the gate's own.
```{r apply-gate}
flips <- unique(impact$compound[impact$sign_changed])
qc_units <- cr_qc_report(ds, gate)
reviewed <- qc_units$well_id[qc_units$verdict == "fail" &
!qc_units$compound %in% flips]
c(gate_excluded = length(gate$excluded), reviewed = length(reviewed))
ds <- cr_apply_gate(ds, gate, units = reviewed)
cr_qc_summary(ds)[, c("step", "cells_before", "cells_after", "percent_removed")]
```
```{r gate-plot}
cr_plot_qc_gate(gate,
statistic = "unit_statistic",
reference = "control_reference",
title = "Each unit against the control of its own batch")
```
## 6. Effect sizes at both levels
The same grid is computed twice from the same cells: once aggregated to the
unit of replication, once on the cells as they are.
```{r effects}
eff_unit <- cr_effect_grid(ds,
value = "log2_fc",
group_var = "treatment",
reference_level = "Vehicle",
by = "compound",
unit = "well_id",
min_n = 3)
eff_cell <- cr_effect_grid(ds,
value = "log2_fc",
group_var = "treatment",
reference_level = "Vehicle",
by = "compound",
min_n = 10)
# shorten the labels of both grids together: cr_compare_levels() joins
# them on the columns they share, so they have to stay consistent
eff_unit$contrast <- sub("^Vehicle -> ", "", eff_unit$contrast)
eff_cell$contrast <- sub("^Vehicle -> ", "", eff_cell$contrast)
eff_unit[, c("compound", "contrast", "n_cmp", "cohens_d",
"cohens_d_ci_low", "cohens_d_ci_high", "p_BH")]
```
`cr_compare_levels()` puts the two side by side and divides the interval
widths.
```{r compare-levels}
cmp <- cr_compare_levels(eff_unit, eff_cell)
head(cmp, 4)
round(attr(cmp, "median_ratio"), 1)
```
The unit-level interval is several times wider than the cell-level interval
for the same estimate. That ratio is not precision gained by measuring more
cells; it is pseudo-replication. Cells within a unit are not independent, so
the cell-level interval understates the uncertainty of the very same number.
Interpretation belongs to the unit level.
Plates are crossed with arms in this design, so the contrast can also be taken
within plate:
```{r blocked}
units <- cr_summarize_wells(ds, channel = "log2_fc", fun = stats::median)
blocked <- cr_blocked_effect(units,
value = "value",
group_var = "treatment",
reference_level = "Vehicle",
block_var = "plate",
by = "compound")
head(blocked[, c("compound", "contrast", "n_blocks", "shift_within_block",
"ci_low", "ci_high", "d_within_block")], 4)
```
```{r screen-figure, fig.height = 4.6}
cr_plot_screen(ds,
value = "log2_fc",
group_var = "treatment",
units = "well_id",
facet_by = "compound",
seed = 1)
```
```{r forest}
cr_plot_forest(eff_unit,
estimate = "cohens_d",
ci_low = "cohens_d_ci_low",
ci_high = "cohens_d_ci_high",
label = "contrast",
facet_by = "compound")
```
## 7. Sample size for the follow-up study
A follow-up powered on the observed point estimate is powered on the luckiest
reading of the current data. `cr_power_grid()` solves the sample size from the
confidence bound *nearer the null* instead, and says so in the `basis` column.
```{r power}
sizes <- cr_power_grid(eff_unit)
sizes[, c("compound", "contrast", "cohens_d", "d_conservative",
"n_observed", "n_conservative")]
```
A contrast whose interval spans the null is not sizable at all, and is
reported as such rather than given a number. The `basis` column says which
of the two happened for every row:
```{r power-basis}
table(sizes$basis)
```
```{r size-figure}
cr_plot_sample_size(sizes, label = "contrast", colour_by = "compound")
```
## 8. The report
`cr_report()` collects the experiment, the quality-control record, the effect
grid, the sample sizes, any supplementary tables and any figures into one
object. Assembling first and formatting later is what keeps a number in the
text tied to the object it came from.
```{r report}
report <- cr_report(
ds,
qc = gate,
effects = eff_unit,
sizes = sizes,
tables = list(units = units, gate_impact = impact, batch_reference = ref),
plots = list(
forest = cr_plot_forest(eff_unit,
estimate = "cohens_d",
ci_low = "cohens_d_ci_low",
ci_high = "cohens_d_ci_high",
label = "contrast",
facet_by = "compound")
),
title = "Four-compound screen",
author = "cellreportR example"
)
names(cr_tables(report))
```
Its `summary` slot is one row per contrast, with the effect grid and the
sample sizes already joined — the overview table of the write-up, derived
rather than retyped:
```{r report-summary}
report$summary[, c("compound", "contrast", "cohens_d", "cohens_d_ci_low",
"cohens_d_ci_high", "n_conservative")]
```
Tables go out as machine-readable files:
```{r export}
out <- file.path(tempdir(), "cr_vignette_out")
dir.create(out, showWarnings = FALSE)
written <- cr_export_tables(cr_tables(report), out, format = "csv",
one_file = FALSE)
basename(written)
```
And every number quoted in prose is emitted from the objects themselves,
rather than transcribed:
```{r macros}
macro_file <- file.path(out, "generated-numbers.tex")
cr_macros(
list(
n_units = nrow(units),
n_cells = cr_n_cells(ds),
top_effect = max(eff_unit$cohens_d),
n_conservative = as.integer(max(sizes$n_conservative, na.rm = TRUE))
),
macro_file
)
cat(readLines(macro_file), sep = "\n")
```
Rendering the assembled object to a document is one further call. It is shown
but not executed here, because it needs a working pandoc and writes a file:
```{r render, eval = FALSE}
cr_render_report(report, output_dir = out, format = "html")
# The same assembled report can also be exported as PDF or Word:
cr_render_report(report, output_dir = out, format = "pdf")
cr_render_report(report, output_dir = out, format = "docx")
```
## What was decided, and where it is recorded
| Decision | Made by | Recorded in |
|---|---|---|
| Which files are one unit | `cr_assign_units()` | `cr_unit_map()`, `provenance` |
| Which objects are debris | `cr_exclude_small()` | `qc_log` |
| How many cells per unit | `cr_balance_cells()` | `qc_log` |
| What a batch is | `batch_vars` | `cr_batch_reference()` |
| Which units are excluded | reviewed list | `qc_log`, `metadata$qc_gate` |
| What the exclusions cost | `cr_qc_gate_impact()` | report table |
| Which level is interpreted | unit | `cr_compare_levels()` |
| What powers the next study | conservative bound | `basis` column |
## Next steps
* `vignette("getting-started")` — the short path on a single plate.
* `vignette("statistical-analysis")` — cell versus unit inference in detail.
* `vignette("dose-response")` — fitting exposure-response curves.
* `vignette("shiny-app")` — the same pipeline in a browser.
