When the spatial data stops fitting
The mif object described in the introduction holds every sample’s spatial data in
memory, as a list of data frames. For TMA cores that is the right thing
to do. A core is a couple of thousand cells, the five cores shipped with
this package come to 12,782 cells in total, and the whole object is a
few megabytes.
Whole slide images are a different problem. A single slide can be a
million cells, and a cohort of them does not fit. The numbers from the
cohort this feature was built for: 283 slides, 342,267,952 cells, 21
columns per cell. As R data frames that is 34.2 GB. On
a 38.7 GB machine the mif is 89% of memory before any
computation starts, which means ripleys_k(workers = 12)
cannot fork a single worker — the run fails before it does any work.
The thing that makes this fixable is that the metrics do not need
most of the data most of the time. ripleys_k() reads a
sample id, four bounding-box columns and however many marker columns you
asked for. On a 21-column HALO export that is 4 columns of 21, and on a
real export with 51 columns it is worse. The other 17 are carried in
memory, in the parent process, for the entire run, and never read.
So the spatial slot can be a manifest of files instead of a list of data frames. Each worker then reads only the columns it needs, for the one sample it is handling, and the parent holds an index rather than a cohort. On those 283 slides the same data is 2.19 GB of parquet against 34.2 GB in memory.
This is not a speed-up
Worth saying plainly, because it is easy to assume otherwise. Disk-backing does not make the metrics faster, and on a single small sample it is slightly slower — there is a file read that was not there before.
What it changes is which memory cost you pay. For Ripley’s K
the per-worker close pair list from the engine overtakes the spatial
frame at around max(r_range) = 30, and at the default
r_range = seq(0, 100, 1) it is roughly twelve times the
frame. The frame is not what makes a large run fail. The parent’s copy
of the whole cohort is, and that is a fixed cost you cannot trade
against workers. Disk-backing removes it; big
is the argument for bounding the pair list.
If your data fits in memory, keep it in memory.
Writing and reading a store
mif_to_disk() takes an ordinary mif and a
directory to write it to, and returns a mif that reads from
that directory.
x <- create_mif(clinical_data = example_clinical %>%
mutate(deidentified_id = as.character(deidentified_id)),
sample_data = example_summary %>%
mutate(deidentified_id = as.character(deidentified_id)),
spatial_list = example_spatial,
patient_id = "deidentified_id",
sample_id = "deidentified_sample")
store <- file.path(store_dir, "cohort.mif")
xd <- mif_to_disk(x, store)
xd
#> 229 patients spanning 229 samples and 5 spatial data frames were foundpath has no default. A store can be tens of gigabytes,
so it is always written where you say it should be and never somewhere
chosen for you — the same reasoning as sigma in
split_tissue().
The result is still a mif, with the same six slots, and
the spatial slot still answers to length(),
names() and [[:
class(xd)
#> [1] "mif"
names(xd)
#> [1] "clinical" "sample" "spatial" "derived" "patient_id"
#> [6] "sample_id"
length(xd$spatial)
#> [1] 5
names(xd$spatial)
#> [1] "TMA1_[3,B].tif" "TMA2_[3,B].tif" "TMA3_[7,B].tif" "TMA3_[9,K].tif"
#> [5] "TMA3_[8,U].tif"
dim(xd$spatial[["TMA3_[9,K].tif"]])
#> [1] 1803 51That last line reads one file off disk. So does $, which
is the convenient way to pull one sample up at the console:
dim(xd$spatial$`TMA3_[9,K].tif`)
#> [1] 1803 51Tab completion works on both, and neither form reads more than it has
to. xd$spatial$Tab offers the sample names and
reads nothing at all — they come from the index.
xd$spatial[["TMA3_[9,K].tif"]]$Tab offers that
sample’s columns, having read that one sample.
Iteration works too, so code written against an in-memory
mif does not need rewriting:
sapply(xd$spatial, nrow)
#> TMA1_[3,B].tif TMA2_[3,B].tif TMA3_[7,B].tif TMA3_[9,K].tif TMA3_[8,U].tif
#> 3803 3008 1850 1803 2318There is a catch worth knowing before you use that on a real cohort.
lapply() and everything built on it materialise their input
before applying anything, so iterating the slot reads
every sample — which is the cost disk-backing exists to
avoid. Iterating over names() instead keeps one sample
resident at a time:
vapply(names(xd$spatial), function(s) nrow(xd$spatial[[s]]), numeric(1))
#> TMA1_[3,B].tif TMA2_[3,B].tif TMA3_[7,B].tif TMA3_[9,K].tif TMA3_[8,U].tif
#> 3803 3008 1850 1803 2318Both give the same answer here, on five cores. On 283 slides the first reads 34 GB and the second reads one slide at a time.
The slot itself is an index:
xd$spatial
#> <mif_store> 5 samples, 12,782 cells
#> root: /var/folders/63/qj4bjln95bs911vg920503rh001s8k/T//RtmpMa6X14/spatialTIME-vignette/cohort.mif
#> TMA1_[3,B].tif (3,803), TMA2_[3,B].tif (3,008), TMA3_[7,B].tif (1,850), ... and 2 more
format(object.size(xd$spatial), units = "KB")
#> [1] "7.9 Kb"
format(object.size(x$spatial), units = "MB")
#> [1] "3.6 Mb"The store is a plain directory you can inspect, copy, archive and move:
list.files(store, recursive = FALSE)
#> [1] "clinical.rds" "derived" "manifest.json" "sample.rds"
#> [5] "spatial"
list.files(file.path(store, "spatial"))
#> [1] "TMA1__3_B_.tif.parquet" "TMA2__3_B_.tif.parquet" "TMA3__7_B_.tif.parquet"
#> [4] "TMA3__8_U_.tif.parquet" "TMA3__9_K_.tif.parquet"manifest.json records the format version, the id
columns, and each sample’s file, row count, byte size and column schema.
Paths inside it are relative, so moving or copying the whole directory
does not break it.
open_mif() reads one back:
Opening validates before returning. Every file is re-probed and compared against the manifest — existence, byte size, row count and column names — so a truncated, half-written or partially copied store fails immediately, with a message naming the sample, rather than somewhere in the middle of a long computation. It costs about a millisecond per sample, because a parquet footer read does not touch the data.
collect_mif() is the inverse of
mif_to_disk(), named after dplyr::collect()
for the same reason: it is where lazily-referenced data gets
materialised.
back <- collect_mif(xd)
class(back$spatial)
#> [1] "list"
all.equal(back$spatial, x$spatial, check.attributes = FALSE)
#> [1] TRUEOn a cohort that does not fit in memory, collecting all of it defeats
the point. samples takes a subset, which is how you pull a
few samples out to look at them:
two <- collect_mif(xd, samples = c("TMA3_[9,K].tif", "TMA1_[3,B].tif"))
names(two$spatial)
#> [1] "TMA3_[9,K].tif" "TMA1_[3,B].tif"
sapply(two$spatial, nrow)
#> TMA3_[9,K].tif TMA1_[3,B].tif
#> 1803 3803Running the metrics
Nothing changes. Every function in the package takes a disk-backed
mif the same way it takes an in-memory one.
mnames <- c("CD3..Opal.570..Positive", "CD8..Opal.520..Positive")
xd <- ripleys_k(mif = xd, mnames = mnames, r_range = seq(0, 100, 10),
edge_correction = "translation", permute = FALSE,
workers = 1, overwrite = TRUE)
k <- xd$derived$univariate_Count
head(k[k$deidentified_sample == "TMA3_[9,K].tif" &
k$Marker == mnames[1],
c("deidentified_sample", "Marker", "r", "Observed K", "Exact CSR")])
#> deidentified_sample Marker r Observed K Exact CSR
#> 67 TMA3_[9,K].tif CD3..Opal.570..Positive 0 0.00000 0.00000
#> 68 TMA3_[9,K].tif CD3..Opal.570..Positive 10 38.53943 18.72003
#> 69 TMA3_[9,K].tif CD3..Opal.570..Positive 20 4651.50907 1268.35887
#> 70 TMA3_[9,K].tif CD3..Opal.570..Positive 30 11112.73225 3188.16333
#> 71 TMA3_[9,K].tif CD3..Opal.570..Positive 40 20118.32327 5871.05676
#> 72 TMA3_[9,K].tif CD3..Opal.570..Positive 50 30962.47935 9204.24046TMA3_[9,K].tif is picked deliberately. The marker counts
in these example cores vary by a factor of thirty — this one has 536
CD3+ cells where TMA1_[3,B].tif has 17 — and on the sparse
cores no two positive cells are within 40 units of each other, so
Observed K is a legitimate but unenlightening column of
zeros.
Derived results are written to the store as they are computed, so they survive a reopen rather than living only in the session that produced them:
z <- open_mif(store)
names(z$derived)
#> [1] "univariate_Count"
identical(z$derived$univariate_Count, xd$derived$univariate_Count)
#> [1] TRUEThe results are the same numbers as the in-memory mif
gives. Not approximately — the same:
x <- ripleys_k(mif = x, mnames = mnames, r_range = seq(0, 100, 10),
edge_correction = "translation", permute = FALSE,
workers = 1, overwrite = TRUE)
all.equal(xd$derived$univariate_Count, x$derived$univariate_Count)
#> [1] TRUEThat is a tested contract rather than an aspiration. The in-memory path ignores the column projection entirely, which means a projection that forgot a column a metric actually uses would fail only on the disk path — so the package carries a parity test comparing the two for every metric, at one and two workers, with and without permutation.
Permuted results match too. Seeds are drawn once in the parent
process and handed to the workers, so a run depends only on your
set.seed() and not on workers:
set.seed(42)
a <- ripleys_k(xd, mnames = mnames[1], r_range = seq(0, 50, 10), permute = TRUE,
num_permutations = 10, workers = 1, overwrite = TRUE)
set.seed(42)
b <- ripleys_k(x, mnames = mnames[1], r_range = seq(0, 50, 10), permute = TRUE,
num_permutations = 10, workers = 1, overwrite = TRUE)
all.equal(a$derived$univariate_Count, b$derived$univariate_Count)
#> [1] TRUEThe store is shared state
That last chunk did something worth noticing. The result went to
a rather than back to xd, but the run still
wrote its results into the store — so the store and the in-memory
xd no longer agree, while the store and a
do:
identical(open_mif(store)$derived$univariate_Count, xd$derived$univariate_Count)
#> [1] FALSE
identical(open_mif(store)$derived$univariate_Count, a$derived$univariate_Count)
#> [1] TRUEThis is the one place a disk-backed mif does not behave
like an ordinary R object. a <- ripleys_k(xd, ...)
leaves the in-memory xd alone, exactly as copy-on-modify
would lead you to expect — but a and xd point
at the same directory, and the metric wrote there. The store reflects
whichever mif most recently ran a metric against it.
Nothing is silently dropped: each run carries forward the derived slots it inherited, so only the slot being written is replaced. But if you want two parameterisations side by side, give them separate stores rather than separate variables:
iso <- mif_to_disk(collect_mif(xd), file.path(store_dir, "iso.mif"))
iso <- ripleys_k(iso, mnames = mnames, r_range = seq(0, 100, 10),
edge_correction = "isotropic", permute = FALSE, overwrite = TRUE)For the ordinary case, assigning the result back to the same name —
xd <- ripleys_k(xd, ...) — keeps the object and the
store in step and the question does not arise.
Pointing at files you already have
If your data is already parquet, create_mif() accepts a
named character vector of file paths instead of a list of data frames.
Nothing is copied.
files <- list.files(file.path(store, "spatial"), full.names = TRUE)
names(files) <- names(xd$spatial)
ref <- create_mif(clinical_data = example_clinical %>%
mutate(deidentified_id = as.character(deidentified_id)),
sample_data = example_summary %>%
mutate(deidentified_id = as.character(deidentified_id)),
spatial_list = files,
patient_id = "deidentified_id",
sample_id = "deidentified_sample")
ref$spatial
#> <mif_store> 5 samples, 12,782 cells
#> referenced files (no store root)
#> TMA1_[3,B].tif (3,803), TMA2_[3,B].tif (3,008), TMA3_[7,B].tif (1,850), ... and 2 moreThis is reference mode, and the print-out says so: there is
no store root, because the files are yours and live wherever you put
them. Each file must hold exactly one sample — every metric uses the
convex hull of all of a file’s cells as the observation window, so a
file with two samples in it would give both the wrong window. If the
vector is unnamed, the names are read from the sample_id
column and a file holding more than one distinct value is rejected.
The trade-off is that nothing may write to those files. A function
that adds a per-cell column has no business modifying your primary data,
so split_tissue() refuses on a reference-mode
mif and tells you to make a store first.
Subsetting
subset_mif() needs somewhere to put the result, for the
same reason mif_to_disk() does: the subset of a cohort too
large for memory usually is too.
sub <- subset_mif(xd, classifier = "Classifier.Label", level = "Tumor",
markers = mnames, path = file.path(store_dir, "tumor.mif"))
cells <- function(m) vapply(names(m$spatial),
function(s) nrow(m$spatial[[s]]), numeric(1))
data.frame(all = cells(xd), tumor_only = cells(sub))
#> all tumor_only
#> TMA1_[3,B].tif 3803 3611
#> TMA2_[3,B].tif 3008 2703
#> TMA3_[7,B].tif 1850 1570
#> TMA3_[9,K].tif 1803 817
#> TMA3_[8,U].tif 2318 1689It writes sample by sample and releases each one, so peak memory stays at a single sample regardless of cohort size.
Checking that the files are still what you think
open_mif() validates at open time, but two situations
are not covered by that. A reference-mode mif has no
manifest, so there was never anything to validate. And a long session
can hold a store open for hours while the files underneath it are
replaced or regenerated.
Both leave the recorded row count stale, and that matters more than
it sounds: every mask in this package is positional, so a stale row
count means markers attached to the wrong cells and plausible, wrong
numbers rather than an error. verify_mif() re-probes
everything and reports the first file that no longer matches.
It returns TRUE invisibly and errors otherwise, so it is
safe to put at the top of a script — wrap it in isTRUE() if
you want to see the answer.
isTRUE(verify_mif(xd))
#> [1] TRUE
isTRUE(verify_mif(ref))
#> [1] TRUE
isTRUE(verify_mif(x)) # in-memory: nothing to check, passes trivially
#> [1] TRUEThat last line is the point of it passing trivially in memory: you do not have to branch on the representation to call it.
Adding columns
split_tissue() writes three per-cell columns, and on a
store those go to a per-sample overlay file rather than into
the base parquet.
xd <- split_tissue(xd, classifier = "Classifier.Label",
class1 = "Tumor", class2 = "Stroma",
sigma = 40, interface_width = 100, workers = 1)
list.files(file.path(store, "overlay"))
#> [1] "TMA1__3_B_.tif.parquet" "TMA2__3_B_.tif.parquet" "TMA3__7_B_.tif.parquet"
#> [4] "TMA3__8_U_.tif.parquet" "TMA3__9_K_.tif.parquet"
# the three columns the overlay added, appended after the base ones
spat <- xd$spatial[["TMA3_[9,K].tif"]]
tail(colnames(spat), 3)
#> [1] "density_compartment" "refined_density_compartment"
#> [3] "density_score"
table(spat$refined_density_compartment)
#>
#> Tumor Interface Stroma
#> 466 882 455Appending columns by rewriting would touch all 2.19 GB of a cohort to
add about 12 MB of new data, so the base files are never rewritten.
Reads cbind the overlay back on, and overlay columns shadow base columns
of the same name, which is what lets a re-run with
overwrite = TRUE replace its own columns instead of
duplicating them.
What survives a round trip, and what does not
Column names, column order, types, and factor levels — including
unused levels and NAs — all round-trip exactly. That last
one matters because split_tissue() writes two factor
columns per cell.
Row names do not. They come back as
1:nrow. A columnar file has nowhere to put them, storing
them would mean an extra column per sample for something no function
here reads, and dplyr resets them anyway. Non-default row
names only arise from having subset a data frame, so this is visible
only if you compare with identical() after doing that. Use
all.equal(check.attributes = FALSE), as above, or reset
them yourself.
A few column types are refused outright rather than stored badly, and the check runs before anything is written:
bad <- example_spatial[["TMA3_[9,K].tif"]][1:5, ]
bad$notes <- as.raw(1:5)
mif_to_disk(toy <- create_mif(
clinical_data = data.frame(deidentified_id = "p1"),
sample_data = data.frame(deidentified_id = "p1",
deidentified_sample = "TMA3_[9,K].tif"),
spatial_list = list(`TMA3_[9,K].tif` = bad),
patient_id = "deidentified_id", sample_id = "deidentified_sample"),
file.path(store_dir, "will-not-exist.mif"))
#> Error: Sample "TMA3_[9,K].tif" has a column that cannot be stored on disk: "notes" (raw).
#> Parquet would not return it unchanged. Drop the column, or keep this mif in memory.raw is the one that justifies the check: parquet accepts
it and hands back an integer, silently. List columns come
back as something that is not what you put in. If you have a column like
that, drop it before writing or keep the mif in memory.
The spatial slot of a store is read-only. The assignment forms that look like they should work on a named list will not, because underneath the slot is a vector of file paths and assigning to it would quietly destroy the index:
xd$spatial[[1]] <- data.frame(a = 1)
#> Error: A disk-backed spatial slot is read-only.
#> Use `collect_mif(mif)` to bring it into memory first, or `split_tissue()` to add per-cell columns through the overlay.Use collect_mif() to bring the data into memory if you
need to modify it directly.
One thing to avoid
Do not resize arrow’s thread pool from inside a worker.
arrow::set_cpu_count() or
set_io_thread_count() called in an mclapply()
child deadlocks — the run hangs at 0% CPU with no error and no output,
which is about the worst failure mode available. Reading parquet in a
child is fine; it is only the resize that hangs.
If you need to limit arrow’s threads, set the count in the parent before calling a metric and the workers will inherit it safely.
arrow::set_cpu_count(2) # in the parent, before the metric
xd <- ripleys_k(xd, mnames = mnames, r_range = 0:100, workers = 6)