TileDB example

.Rbuildignore

@@ -16,3 +16,4 @@
 ^Apptainer\.def$
 ^README\.Rmd$
 ^README\.html$
+^binder$

.lintr

@@ -8,6 +8,7 @@ linters: linters_with_defaults(
   )
 exclusions: list(
     "inst/doc",
+    "inst/examples",
     "man",
     "vignettes",
     "data-raw",

DESCRIPTION

@@ -50,7 +50,9 @@ Suggests:
     lintr,
     stringr,
     tiledb,
-    TileDBArray
+    TileDBArray,
+    cowplot,
+    png
 Config/testthat/edition: 3
 VignetteBuilder: knitr
 URL: https://pennlinc.github.io/ModelArray, https://pennlinc.github.io/ModelArray/

R/utils.R

@@ -73,6 +73,11 @@ flagAnalysisExistInh5 <- function(fn_h5, analysis_name) {
   if (grepl("\\.tdb/?$", filepath, ignore.case = TRUE)) {
     return("tiledb")
   }
+  # A remote URI (s3://, gcs://, azure://, ...) can only be a TileDB store;
+  # rhdf5 reads local files only. file:// is left to the local checks below.
+  if (.is_remote_uri(filepath)) {
+    return("tiledb")
+  }
   # Otherwise only claim a directory when it actually holds a store. Writers
   # delete the group they are about to replace, so a mistyped output directory
   # must not be mistaken for one.
@@ -269,6 +274,45 @@ flagAnalysisExistInh5 <- function(fn_h5, analysis_name) {
 }
 
 
+#' Is this a URI with a scheme other than file:// (for example s3://)?
+#' @noRd
+.is_remote_uri <- function(path) {
+  grepl("^[[:alpha:]][[:alnum:]+.-]*://", path) && !grepl("^file://", path, ignore.case = TRUE)
+}
+
+
+#' A TileDB VFS handle on the current default context
+#'
+#' tiledb::tiledb_get_vfs() caches a VFS built from the first context it saw,
+#' which would miss S3 settings a user applied later with tiledb_ctx().
+#' @noRd
+.tiledb_vfs <- function() {
+  tiledb::tiledb_vfs(ctx = tiledb::tiledb_get_context())
+}
+
+
+#' Does a TileDB store root exist, as a local directory or as a URI?
+#' @noRd
+.tiledb_store_exists <- function(uri) {
+  if (dir.exists(uri)) {
+    return(TRUE)
+  }
+  isTRUE(tryCatch(tiledb::tiledb_vfs_is_dir(uri, vfs = .tiledb_vfs()), error = function(e) FALSE))
+}
+
+
+#' Names of the entries directly under a TileDB URI (character(0) when absent)
+#'
+#' Uses TileDB's VFS, so it works for local paths, file:// and remote URIs
+#' alike. Plain directories are listed too (test fixtures store arrays under
+#' them), so callers still check for the array they need.
+#' @noRd
+.tiledb_list_children <- function(uri) {
+  children <- tryCatch(tiledb::tiledb_vfs_ls(uri, vfs = .tiledb_vfs()), error = function(e) character(0))
+  sort(unique(basename(sub("/+$", "", children))))
+}
+
+
 #' Convert TileDB metadata values to a character vector
 #' @noRd
 .parse_tiledb_names <- function(x) {
@@ -1126,12 +1170,28 @@ bind_cols_check_emptyTibble <- function(a, b) {
 #' into memory. The returned
 #' object has a \code{print} method that displays a formatted summary.
 #'
+#' A TileDB store can also be remote, for example in an S3 bucket. Configure
+#' the TileDB context before calling, exactly as you would before
+#' \code{\link{ModelArray}()}:
+#' \preformatted{
+#' cfg <- tiledb::tiledb_config()
+#' cfg["vfs.s3.region"] <- "us-east-1"
+#' cfg["vfs.s3.no_sign_request"] <- "true"   # public bucket, no credentials
+#' tiledb::tiledb_ctx(cfg)
+#' ModelArraySummary("s3://bucket/path/store.tdb")
+#' }
+#' The same settings can come from environment variables
+#' (\code{TILEDB_VFS_S3_REGION}, \code{TILEDB_VFS_S3_NO_SIGN_REQUEST}).
+#' Each scalar costs a few requests, so summarizing a store with many scalars
+#' over the network takes seconds rather than milliseconds.
+#'
 #' \code{h5summary()} is a backward-compatible alias for
 #' \code{ModelArraySummary()}. Despite the historical name, it supports both
 #' HDF5 files and TileDB stores.
 #'
 #' @param filepath Character. Path to an HDF5 (\code{.h5}) file or TileDB
-#'   (\code{.tdb}) store.
+#'   (\code{.tdb}) store. A TileDB store may be given as a URI
+#'   (\code{s3://...}, \code{file://...}); HDF5 files must be local paths.
 #' @param backend Character. Storage backend: \code{"auto"} (default),
 #'   \code{"hdf5"}, or \code{"tiledb"}.
 #'
@@ -1155,6 +1215,7 @@ bind_cols_check_emptyTibble <- function(a, b) {
 #' \dontrun{
 #' ModelArraySummary("path/to/data.h5")
 #' ModelArraySummary("path/to/store.tdb")
+#' ModelArraySummary("s3://bucket/path/store.tdb") # after configuring the TileDB context
 #'
 #' # Inspect before deciding which scalars to load
 #' info <- ModelArraySummary("path/to/data.h5")
@@ -1185,18 +1246,15 @@ h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
 
 #' @noRd
 .modelarray_summary_impl <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
-  if (!file.exists(filepath)) {
-    stop("File not found: ", filepath)
-  }
   backend <- .resolve_storage_backend(filepath, backend)
   if (identical(backend, "tiledb")) {
     .require_tiledb_support()
-    scalars_dir <- .modelarray_path(filepath, "scalars")
-    scalar_candidates <- if (dir.exists(scalars_dir)) {
-      basename(list.dirs(scalars_dir, recursive = FALSE, full.names = TRUE))
-    } else {
-      character(0)
+    # TileDB stores may live behind a URI (s3://, file://, ...), so existence and
+    # listing go through TileDB's VFS rather than file.exists()/list.dirs().
+    if (!.tiledb_store_exists(filepath)) {
+      stop("TileDB store not found: ", filepath)
     }
+    scalar_candidates <- .tiledb_list_children(.modelarray_path(filepath, "scalars"))
     scalar_names <- scalar_candidates[vapply(
       scalar_candidates,
       function(scalar_name) {
@@ -1240,12 +1298,7 @@ h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
       )
     }
 
-    results_dir <- .modelarray_path(filepath, "results")
-    analysis_candidates <- if (dir.exists(results_dir)) {
-      basename(list.dirs(results_dir, recursive = FALSE, full.names = TRUE))
-    } else {
-      character(0)
-    }
+    analysis_candidates <- .tiledb_list_children(.modelarray_path(filepath, "results"))
     analyses <- analysis_candidates[vapply(
       analysis_candidates,
       function(analysis_name) {
@@ -1267,6 +1320,9 @@ h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
     )
   }
 
+  if (!file.exists(filepath)) {
+    stop("File not found: ", filepath)
+  }
   listing <- rhdf5::h5ls(filepath)
   rhdf5::h5closeAll()
 

README.Rmd

@@ -49,6 +49,7 @@ That's it.
 **It scales.** ModelArray uses [HDF5](https://www.hdfgroup.org/solutions/hdf5/) by default, and can also read TileDB stores when the optional [TileDB](https://tiledb.com/) R packages are installed. Both backends use [DelayedArray](https://bioconductor.org/packages/DelayedArray/) for lazy access.
 An ABCC dMRI dataset has ~350,000 voxels × ~26,000-sessions - about 291 GB if loaded in its entirety - never enters RAM.
 ModelArray reads one element at a time, so memory usage stays flat regardless of dataset size.
+TileDB stores can also live in the cloud: `vignette("rbc-s3-tiledb")` fits GAMs to the 553 GB RBC cortical thickness store straight from a public S3 bucket, and you can run it in your browser with [Binder](https://mybinder.org/v2/gh/PennLINC/ModelArray/tiledb-example?urlpath=rstudio).
 
 **It's flexible.** Beyond linear models, ModelArray supports:
 

README.md

@@ -50,6 +50,10 @@ access. An ABCC dMRI dataset has ~350,000 voxels × ~26,000-sessions -
 about 291 GB if loaded in its entirety - never enters RAM. ModelArray
 reads one element at a time, so memory usage stays flat regardless of
 dataset size.
+TileDB stores can also live in the cloud: `vignette("rbc-s3-tiledb")`
+fits GAMs to the 553 GB RBC cortical thickness store straight from a
+public S3 bucket, and you can run it in your browser with
+[Binder](https://mybinder.org/v2/gh/PennLINC/ModelArray/tiledb-example?urlpath=rstudio).
 
 **It’s flexible.** Beyond linear models, ModelArray supports:
 

_pkgdown.yml

@@ -34,6 +34,8 @@ navbar:
           href: articles/modelling.html
         - text: "Exploring Data with Convenience Functions"
           href: articles/exploring-modelarray-data.html
+        - text: "Cloud TileDB: RBC Cortical Thickness from S3"
+          href: articles/rbc-s3-tiledb.html
         - text: "-------"
         - text: Resources
         - text: "Containers"

binder/apt.txt

@@ -0,0 +1,6 @@
+libhdf5-dev
+libcurl4-openssl-dev
+libssl-dev
+libxml2-dev
+libgit2-dev
+bzip2

binder/install.R

@@ -0,0 +1,10 @@
+# Binder image for vignette("rbc-s3-tiledb"): everything the example report needs.
+# CRAN packages come from the Posit Package Manager snapshot named in runtime.txt.
+install.packages(c("remotes", "BiocManager", "tiledb", "mgcv", "ggplot2", "gratia",
+                   "gamm4", "lme4", "MASS", "rmarkdown", "knitr", "jsonlite", "cowplot", "png"))
+BiocManager::install(c("TileDBArray", "HDF5Array", "rhdf5", "DelayedArray"),
+                     update = FALSE, ask = FALSE)
+options(repos = BiocManager::repositories())
+remotes::install_github("andy1764/ComBatFamily", upgrade = "never")
+# This repository is ModelArray itself, at the branch Binder was launched from.
+remotes::install_local(".", upgrade = "never", force = TRUE)

binder/postBuild

@@ -0,0 +1,13 @@
+#!/usr/bin/env bash
+# Runs once when the Binder image is built (as the notebook user, no root).
+set -euo pipefail
+
+# Connectome Workbench from conda-forge, installed into the home directory with micromamba.
+mkdir -p "$HOME/.local/bin"
+curl -Ls https://micro.mamba.pm/api/micromamba/linux-64/latest | tar -xj -C "$HOME/.local" bin/micromamba
+"$HOME/.local/bin/micromamba" create -y -p "$HOME/.wb" -c conda-forge connectome-workbench
+"$HOME/.local/bin/micromamba" clean -y --all
+
+# Put the report where the RStudio session opens, so "Knit" is the only click needed.
+mkdir -p "$HOME/rbc_s3_example"
+cp inst/examples/rbc-s3-tiledb/rbc_s3_tiledb_example.Rmd "$HOME/rbc_s3_example/"

binder/runtime.txt

@@ -0,0 +1 @@
+r-4.4-2025-03-01

binder/start

@@ -0,0 +1,7 @@
+#!/usr/bin/env bash
+# Runs before the server starts in every Binder session.
+export PATH="$HOME/.wb/bin:$PATH"
+# Anonymous reads from the public RBC bucket; the report also sets these itself.
+export TILEDB_VFS_S3_NO_SIGN_REQUEST=true
+export TILEDB_VFS_S3_REGION=us-east-1
+exec "$@"

inst/examples/rbc-s3-tiledb/README.md

@@ -0,0 +1,38 @@
+# RBC cortical thickness from S3: a runnable ModelArray example
+
+`rbc_s3_tiledb_example.Rmd` reads the RBC FreeSurferPost TileDB store straight from
+its public S3 bucket, fits vertex-wise GAMs of cortical thickness on age in a handful
+of Glasser parcels, and writes the results as CIFTI files you can open in Connectome
+Workbench. `vignette("rbc-s3-tiledb")` explains how to run it and what to expect.
+
+## Files
+
+| File | What it is |
+|---|---|
+| `rbc_s3_tiledb_example.Rmd` | The report. Render it with `rmarkdown::render()`; see the vignette for the exact call and parameters. |
+| `glasser360_fsLR164k_labels.rds` | One Glasser (HCP-MMP1.0) parcel key per column of the TileDB store (327,684 fsLR-164k vertices, left hemisphere then right, medial wall = 0), with the label table and each parcel's sensorimotor-association (S-A) axis rank. 26 KB. |
+| `sa_axis_glasser360_ranks.csv` | S-A axis average ranks for the 180 Glasser regions (`final.rank`: 1 = sensorimotor pole, 180 = association pole). |
+| `make_glasser_fsLR164k_labels.sh`, `make_glasser_fsLR164k_labels.R` | How the `.rds` was made: the 32k Glasser dlabel resampled to 164k with Connectome Workbench and TemplateFlow's fsLR spheres. Kept for provenance; you don't need to run them. |
+| `brain_sa_axis.png`, `brain_sa_3b.png`, `brain_sa_46.png` | Lateral views of the left inflated cortex coloured by the S-A axis in the colour scheme of Sydnor et al. 2021 Figure 2A (gold = sensorimotor, purple = association). `brain_sa_axis` outlines and labels both areas 3b and 46; the other two outline one parcel each. The figures use them as locators. |
+| `make_brain_sa_axis.R` | How the three brain PNGs were made (filled-mesh render of the 164k inflated surface). Provenance; you don't need to run it. |
+
+## Sources and attribution
+
+- **Parcellation:** Glasser MF et al. (2016). A multi-modal parcellation of human
+  cerebral cortex. *Nature* 536, 171–178. The 32k dlabel is the copy distributed with
+  [XCP-D](https://github.com/PennLINC/xcp_d) (`atlas-Glasser_space-fsLR_den-32k_dseg.dlabel.nii`);
+  the HCP-MMP1.0 parcellation is released under the HCP open-access data use terms.
+- **S-A axis ranks:** Sydnor VJ et al. (2021). Neurodevelopment of the association
+  cortices: patterns, mechanisms, and implications for psychopathology. *Neuron* 109,
+  2820–2846. File from [PennLINC/S-A_ArchetypalAxis](https://github.com/PennLINC/S-A_ArchetypalAxis)
+  (`Glasser360_MMP/Sensorimotor_Association_Axis_AverageRanks.csv`). Please cite
+  Sydnor et al. 2021 if you use the ranks. The S-A axis brain (`brain_sa_axis.png`)
+  reuses that paper's Figure 2A colour scheme (gold → white → purple, diverging at
+  rank 90).
+- **Surfaces for resampling and display:** [TemplateFlow](https://www.templateflow.org/)
+  `tpl-fsLR` (spheres, vertex-area maps and 32k display surfaces), downloaded by the
+  scripts when needed.
+- **Data:** Reproducible Brain Charts (RBC) FreeSurferPost release, published as a
+  TileDB store at `s3://fcp-indi/data/Projects/RBC/RBC_FreeSurfer-Post_ModelArray/`.
+  See the RBC documentation at <https://reprobrainchart.github.io/> for the data
+  release and its terms.

inst/examples/rbc-s3-tiledb/make_brain_sa_axis.R

@@ -0,0 +1,95 @@
+#!/usr/bin/env Rscript
+# Build the S-A axis brain images that ship with the example: lateral views of the left
+# inflated fsLR-164k surface coloured by the Glasser sensorimotor-association (S-A) axis,
+# in the colour scheme of Sydnor et al. 2021 (Neuron), Figure 2A -- gold (sensorimotor)
+# through white to purple (association), diverging at rank 90.
+#
+#   make_brain_sa_axis.R <L inflated 164k surf.gii> <glasser360_fsLR164k_labels.rds> <out dir> [wb_command]
+#
+# Writes three PNGs into <out dir>:
+#   brain_sa_axis.png  areas 3b and 46 outlined and labelled (the main locator)
+#   brain_sa_3b.png    area 3b outlined only (for the 3b panel of the vertex figure)
+#   brain_sa_46.png    area 46 outlined only (for the 46 panel of the vertex figure)
+#
+# The 164k inflated surface is the one the example builds from TemplateFlow's 32k surfaces
+# (see the report's cifti display step, or make_glasser_fsLR164k_labels.sh). This is
+# provenance for files that ship with the example; users do not run it. Needs wb_command
+# (to convert the GIFTI to ASCII), xml2, ggplot2 and scales.
+#
+# S-A axis ranks: Sydnor VJ et al. 2021, Neuron (PennLINC/S-A_ArchetypalAxis). Cite Sydnor
+# et al. 2021 if you use them. The surface is rendered as shaded, back-face-culled filled
+# triangles; parcel outlines are the mesh edges between in- and out-of-parcel triangles.
+suppressMessages({library(xml2); library(ggplot2); library(scales)})
+args <- commandArgs(trailingOnly = TRUE)
+if (length(args) < 3) stop("usage: make_brain_sa_axis.R L_inflated.surf.gii labels.rds out_dir [wb_command]")
+surf_gii <- args[1]; labels_rds <- args[2]; out_dir <- args[3]
+wb <- if (length(args) >= 4) args[4] else "wb_command"
+dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
+
+ascii <- tempfile(fileext = ".surf.gii")
+stopifnot(system2(wb, c("-gifti-convert", "ASCII", surf_gii, ascii)) == 0)
+doc <- read_xml(ascii)
+V <- matrix(as.numeric(scan(text = xml_text(xml_find_first(doc, ".//DataArray[@Intent='NIFTI_INTENT_POINTSET']//Data")), quiet = TRUE)), ncol = 3, byrow = TRUE)
+Tt <- matrix(as.integer(scan(text = xml_text(xml_find_first(doc, ".//DataArray[@Intent='NIFTI_INTENT_TRIANGLE']//Data")), quiet = TRUE)), ncol = 3, byrow = TRUE) + 1L
+labels <- readRDS(labels_rds)
+keyL <- labels$key[seq_len(nrow(V))]
+stopifnot(nrow(V) == 163842L)
+
+SA_COLS <- c("goldenrod1", "white", "#6f1282"); SA_MID <- 90; SA_RANGE <- c(1, 180)   # Sydnor 2021 Fig 2A
+rank_col <- function(r) gradient_n_pal(SA_COLS)(rescale_mid(r, to = c(0, 1), from = SA_RANGE, mid = SA_MID))
+WALL <- "grey92"
+rankv <- labels$table$sa_rank[match(keyL, labels$table$key)]
+sacol <- ifelse(is.na(rankv), WALL, rank_col(rankv))
+
+blend <- function(hex, f) { m <- grDevices::col2rgb(hex) / 255; m <- sweep(m, 2, f, `*`); m[m > 1] <- 1; m[m < 0] <- 0; grDevices::rgb(m[1,], m[2,], m[3,]) }
+# left lateral view: anterior (+Y) to the left, superior (+Z) up; keep lateral-facing triangles
+xi <- -V[, 2]; yi <- V[, 3]
+a <- V[Tt[,1],]; b <- V[Tt[,2],]; c <- V[Tt[,3],]
+n <- cbind((b[,2]-a[,2])*(c[,3]-a[,3])-(b[,3]-a[,3])*(c[,2]-a[,2]),
+           (b[,3]-a[,3])*(c[,1]-a[,1])-(b[,1]-a[,1])*(c[,3]-a[,3]),
+           (b[,1]-a[,1])*(c[,2]-a[,2])-(b[,2]-a[,2])*(c[,1]-a[,1]))
+nn <- sqrt(rowSums(n^2)); nn[nn == 0] <- 1; n <- n / nn
+ff <- n[, 1] < 0; Tf <- Tt[ff, , drop = FALSE]; nf <- n[ff, , drop = FALSE]
+light <- c(-0.5, 0.25, 0.85); light <- light / sqrt(sum(light^2))
+shade <- pmax(0.6, 0.78 + 0.22 * as.vector(nf %*% light))
+tcol <- blend(sacol[Tf[, 1]], shade)
+depth <- (V[Tf[,1],1] + V[Tf[,2],1] + V[Tf[,3],1]) / 3
+ord <- order(depth, decreasing = TRUE); Tf <- Tf[ord, ]; tcol <- tcol[ord]; m <- nrow(Tf)
+poly <- data.frame(x = c(xi[Tf[,1]], xi[Tf[,2]], xi[Tf[,3]]), y = c(yi[Tf[,1]], yi[Tf[,2]], yi[Tf[,3]]),
+                   grp = rep(seq_len(m), 3), fill = rep(tcol, 3))
+poly <- poly[order(poly$grp), ]
+base_plot <- ggplot(poly, aes(x, y, group = grp, fill = fill)) + geom_polygon(colour = NA) +
+  scale_fill_identity() + coord_fixed() + theme_void()
+
+outline_df <- function(nm) {
+  inpar <- keyL == labels$table$key[labels$table$name == nm]
+  tri_in <- rowSums(matrix(inpar[Tf], ncol = 3)) >= 2L
+  e <- rbind(Tf[,c(1,2)], Tf[,c(2,3)], Tf[,c(1,3)]); e <- cbind(pmin(e[,1],e[,2]), pmax(e[,1],e[,2]))
+  k <- paste(e[,1], e[,2])
+  keep <- ave(rep(1L, nrow(e)), k, FUN = sum) == 2L & ave(as.integer(rep(tri_in, 3)), k, FUN = sum) == 1L & !duplicated(k)
+  e <- e[keep, , drop = FALSE]
+  data.frame(x = xi[e[,1]], y = yi[e[,1]], xend = xi[e[,2]], yend = yi[e[,2]])
+}
+centroid <- function(nm) { v <- which(keyL == labels$table$key[labels$table$name == nm]); data.frame(x = mean(xi[v]), y = mean(yi[v])) }
+
+# region name -> (glasser parcel, label text, label nudge)
+PARCELS <- list("3b" = list(roi = "L_3b_ROI", dx = -14, dy = 16),
+                "46" = list(roi = "L_46_ROI", dx = -2,  dy = 18))
+render <- function(regions, label) {
+  g <- base_plot
+  for (rg in regions) g <- g + geom_segment(data = outline_df(PARCELS[[rg]]$roi),
+    aes(x, y, xend = xend, yend = yend), inherit.aes = FALSE, colour = "black", linewidth = 0.45)
+  if (label) {
+    labs <- do.call(rbind, lapply(regions, function(rg) cbind(centroid(PARCELS[[rg]]$roi), lab = rg,
+      dx = PARCELS[[rg]]$dx, dy = PARCELS[[rg]]$dy)))
+    g <- g + geom_text(data = labs, aes(x + dx, y + dy, label = lab), inherit.aes = FALSE, colour = "black", fontface = 2, size = 5.4)
+  }
+  g
+}
+save_png <- function(g, file) {
+  ggsave(file.path(out_dir, file), g, width = 3.4, height = 2.8, dpi = 200, bg = "transparent")
+  message(sprintf("wrote %s (%.0f KB)", file, file.size(file.path(out_dir, file)) / 1024))
+}
+save_png(render(c("3b", "46"), label = TRUE),  "brain_sa_axis.png")
+save_png(render("3b", label = FALSE),          "brain_sa_3b.png")
+save_png(render("46", label = FALSE),          "brain_sa_46.png")