From 94d450c9d0a8d7c5aee86af5636dbb5a3e07f975 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 1 Sep 2026 09:03:51 +0930 Subject: [PATCH] add qmd for smooths --- .../index.qmd | 350 ++++++++++++++++++ 1 file changed, 350 insertions(+) create mode 100644 posts/2026-08-06-sccomp-smooth-effects/index.qmd diff --git a/posts/2026-08-06-sccomp-smooth-effects/index.qmd b/posts/2026-08-06-sccomp-smooth-effects/index.qmd new file mode 100644 index 00000000..7381b976 --- /dev/null +++ b/posts/2026-08-06-sccomp-smooth-effects/index.qmd @@ -0,0 +1,350 @@ +--- +title: "sccomp now supports smooth effects" +author: "Stefano Mangiola" +contributors: + - Chen Zhan +date: "2026-08-06" +package: sccomp +tags: + - tidyomics/tidyomicsBlog + - sccomp + - single-cell + - compositional-data + - splines + - pseudotime +description: "Model non-linear changes in cell-type composition with smooth terms in sccomp." +image: cellnexus-plasma-age-smooth.png +format: + html: + toc: true + toc-float: true +execute: + freeze: true + message: false + warning: false +--- + +Continuous biological processes are rarely perfectly linear. Cell-type proportions can rise, plateau, or peak along pseudotime, age, dose, or time since treatment. A straight-line effect can miss these patterns or describe them poorly. + +We are pleased to announce that [`sccomp`](https://mangiolalaboratory.github.io/sccomp/) now supports **smooth terms** for modelling non-linear effects on cell-type composition. + +To demonstrate the new functionality with real data, we use [`cellNexus`](https://github.com/MangiolaLaboratory/cellNexus) metadata to ask how immune-cell composition in blood changes across age. + +*The quality-control workflow was adapted from an analysis by contributor Chen Zhan.* + +The current public `cellNexus` metadata table used for this post does not expose a disease field. Therefore, unlike the source QC notebook, this reproducible example cannot apply `disease == "normal"` and should not be interpreted as a healthy-only analysis. + +# Prepare the CellNexus data + +`cellNexus::get_metadata()` returns a lazy table, so filtering happens before the result is collected into memory. We retain live singlets from blood, remove empty droplets, require age metadata, and keep a defined set of immune cell types. + +```{r} +#| label: prepare-cellnexus-metadata + +library(cellNexus) +library(tidyverse) +library(purrr) +library(sccomp) + +metadata <- get_metadata() |> + filter( + !empty_droplet, + alive, + scDblFinder.class == "singlet", + !is.na(age_days), + tissue_groups == "blood", + age_days > 365 * 10, + cell_type_unified_ensemble %in% c( + "cd8 naive", "cd16 mono", "cd4 tcm", "cd4 th17 em", + "granulocyte", "cd4 th1/th17 em", "treg", "b memory", + "b naive", "nk", "plasma", "cd4 th2 em", "mast", + "cd4 th1 em", "cd8 tem", "mait", "tgd", "cdc", + "cd4 fh em", "cd4 naive", "nkt", "macrophage", + "cd8 tcm", "cd14 mono", "pdc", "ilc" + ) + ) |> + mutate( + age_days_scaled = age_days / sd(age_days), + sample = paste(dataset_id, donor_id, age_days, sep = "__") + ) +``` + +# QC + +```{r} +#| label: qc-cellnexus-metadata + +datasets_to_keep <- metadata |> + count( + dataset_id, + cell_type_unified_ensemble, + name = "n" + ) |> + collect() |> + group_by(dataset_id) |> + mutate(proportion = n / sum(n)) |> + summarise( + entropy = -sum(proportion * log(proportion)), + largest_cell_type_proportion = max(proportion), + .groups = "drop" + ) |> + filter( + entropy > 0, + largest_cell_type_proportion < 0.95 + ) |> + pull(dataset_id) + +metadata <- metadata |> + filter(dataset_id %in% datasets_to_keep) +``` + +The entropy check follows the enrichment QC used in the source analysis: a dataset containing only one annotated cell type has zero entropy. We also exclude datasets in which a single cell type accounts for at least 95% of cells. This is an operational screen for targeted or strongly enriched data, not proof that every retained sample was generated without enrichment. + +We keep all donors that pass this screen, rather than subsampling, so that the blood age trend can use as much information as possible. + +# Fit a five-component smooth + +The new functionality uses the `s()` syntax familiar from [`mgcv`](https://cran.r-project.org/package=mgcv) and [`brms`](https://paulbuerkner.com/brms/). Restricting to blood lets us fit a single shared age smooth without tissue factor smooths. A dataset random intercept still accounts for study-level baseline differences across the many blood datasets. + +```{r} +#| label: fit-age-smooth +#| cache: true +#| eval: false + +fit <- metadata |> + collect() |> + sccomp_estimate( + formula_composition = + ~ s(age_days_scaled, k = 5) + + (1 | dataset_id), + sample = "sample", + cell_group = "cell_type_unified_ensemble", + max_sampling_iterations = 1000, + inference_method = "hmc", + mcmc_seed = 42 + ) + +fit |> saveRDS("fit.rds") +``` + +```{r} +#| label: load-fit + +fit <- readRDS("fit.rds") +``` + +```{r} +#| label: prepare-fit-summary +#| include: false + +model_counts <- attr(fit, "count_data") +``` + +This run included `r scales::comma(n_distinct(model_counts$sample))` blood samples from `r scales::comma(n_distinct(model_counts$dataset_id))` datasets. It modelled `r scales::comma(sum(model_counts$count))` immune cells, including `r scales::comma(sum(model_counts$count[model_counts$cell_type_unified_ensemble == "plasma"]))` plasma cells, `r scales::comma(sum(model_counts$count[model_counts$cell_type_unified_ensemble == "cd4 naive"]))` CD4 naïve cells, and `r scales::comma(sum(model_counts$count[model_counts$cell_type_unified_ensemble == "cd8 naive"]))` CD8 naïve cells. + +# Predict and visualise the curve + +We keep only the age smooth for visualisation: `sccomp_remove_unwanted_effects()` factors out dataset effects from the observed proportions, and `sccomp_predict()` draws the age-only curve over a simple sweep of scaled ages from the fitted data. + +```{r} +#| label: predict-age-smooth + +adjusted <- fit |> + sccomp_remove_unwanted_effects( + formula_composition_keep = ~ s(age_days_scaled, k = 5) + ) |> + left_join( + model_counts |> distinct(sample, age_days_scaled), + by = "sample" + ) + +prediction_grid <- tibble( + age_days_scaled = seq( + min(model_counts$age_days_scaled), + max(model_counts$age_days_scaled), + length.out = 80 + ), + sample = paste0("grid_", 1:80) +) + +curves <- fit |> + sccomp_predict( + formula_composition = ~ s(age_days_scaled, k = 5), + new_data = prediction_grid, + number_of_draws = 500, + summary_instead_of_draws = TRUE + ) |> + arrange(age_days_scaled) + +plasma_adjusted <- adjusted |> + filter(cell_type_unified_ensemble == "plasma") + +plasma_curve <- curves |> + filter(cell_type_unified_ensemble == "plasma") + +cd4_naive_adjusted <- adjusted |> + filter(cell_type_unified_ensemble == "cd4 naive") + +cd4_naive_curve <- curves |> + filter(cell_type_unified_ensemble == "cd4 naive") + +cd8_naive_adjusted <- adjusted |> + filter(cell_type_unified_ensemble == "cd8 naive") + +cd8_naive_curve <- curves |> + filter(cell_type_unified_ensemble == "cd8 naive") +``` + +```{r} +#| label: plot-helpers +#| include: false + +arcsine_sqrt_transform <- scales::new_transform( + name = "arcsine-sqrt", + transform = \(x) asin(sqrt(x)), + inverse = \(x) sin(x)^2, + domain = c(0, 1) +) + +plot_age_smooth <- function(adjusted_data, curve_data, cell_label) { + ggplot() + + geom_point( + data = adjusted_data, + aes(age_days_scaled, adjusted_proportion), + colour = "#5B4B8A", + alpha = 0.25, + size = 1 + ) + + geom_ribbon( + data = curve_data, + aes( + age_days_scaled, + ymin = proportion_lower, + ymax = proportion_upper, + group = 1 + ), + fill = "#D9CBE8", + alpha = 0.75 + ) + + geom_line( + data = curve_data, + aes(age_days_scaled, proportion_mean, group = 1), + colour = "#5B4B8A", + linewidth = 1 + ) + + scale_y_continuous( + transform = arcsine_sqrt_transform, + labels = scales::label_percent() + ) + + labs( + x = "Scaled age", + y = paste0( + cell_label, + " proportion among blood immune cells (arcsine–square-root scale)" + ), + title = paste0( + cell_label, + " composition in blood across age" + ), + subtitle = "sccomp smooth (k = 5); dataset effects removed" + ) + + theme_minimal(base_size = 11) +} +``` + +```{r} +#| label: plot-age-smooth +#| fig-width: 8 +#| fig-height: 5 +#| fig-alt: "Scatter plot of adjusted plasma-cell proportions against scaled age in blood. Points are sample proportions after removing dataset effects. A purple sccomp smooth shows the estimated non-linear age trend, with a light-purple uncertainty ribbon." + +plasma_plot <- plot_age_smooth( + plasma_adjusted, + plasma_curve, + "Plasma-cell" +) + +plasma_plot +``` + +```{r} +#| label: plot-cd4-naive-smooth +#| fig-width: 8 +#| fig-height: 5 +#| fig-alt: "Scatter plot of adjusted CD4 naive proportions against scaled age in blood. Points are sample proportions after removing dataset effects. A purple sccomp smooth shows the estimated non-linear age trend, with a light-purple uncertainty ribbon." + +cd4_naive_plot <- plot_age_smooth( + cd4_naive_adjusted, + cd4_naive_curve, + "CD4 naïve" +) + +cd4_naive_plot +``` + +```{r} +#| label: plot-cd8-naive-smooth +#| fig-width: 8 +#| fig-height: 5 +#| fig-alt: "Scatter plot of adjusted CD8 naive proportions against scaled age in blood. Points are sample proportions after removing dataset effects. A purple sccomp smooth shows the estimated non-linear age trend, with a light-purple uncertainty ribbon." + +cd8_naive_plot <- plot_age_smooth( + cd8_naive_adjusted, + cd8_naive_curve, + "CD8 naïve" +) + +cd8_naive_plot +``` + +```{r} +ggsave( + "cellnexus-plasma-age-smooth.png", + plasma_plot, + width = 8, + height = 5, + dpi = 180 +) + +ggsave( + "cellnexus-cd4-naive-age-smooth.png", + cd4_naive_plot, + width = 8, + height = 5, + dpi = 180 +) + +ggsave( + "cellnexus-cd8-naive-age-smooth.png", + cd8_naive_plot, + width = 8, + height = 5, + dpi = 180 +) +``` + +For plasma cells, the fitted age-only mean starts near `r scales::percent(first(plasma_curve$proportion_mean), accuracy = 0.1)`, reaches a maximum of `r scales::percent(max(plasma_curve$proportion_mean), accuracy = 0.1)`, and ends near `r scales::percent(last(plasma_curve$proportion_mean), accuracy = 0.1)`. For CD4 naïve cells, the corresponding values are `r scales::percent(first(cd4_naive_curve$proportion_mean), accuracy = 0.1)`, `r scales::percent(max(cd4_naive_curve$proportion_mean), accuracy = 0.1)`, and `r scales::percent(last(cd4_naive_curve$proportion_mean), accuracy = 0.1)`. For CD8 naïve cells, they are `r scales::percent(first(cd8_naive_curve$proportion_mean), accuracy = 0.1)`, `r scales::percent(max(cd8_naive_curve$proportion_mean), accuracy = 0.1)`, and `r scales::percent(last(cd8_naive_curve$proportion_mean), accuracy = 0.1)`. Points show sample proportions after removing dataset effects. This is an exploratory demonstration of the smooth interface, not evidence for a general biological age effect. + +# What the smooth is doing + +Under the hood, `sccomp` separates `s(age_days_scaled, k = 5)` into an unpenalised linear component and penalised non-linear basis components. Scaling age improves numerical stability without changing the shape shown on the original year scale. The individual basis coefficients are not usually meaningful on their own; the fitted curve is the quantity to interpret. Decreasing `k` produces a simpler curve, while increasing it gives the model more capacity to bend. + +# Learn more + +The complete [smooth-terms vignette](https://mangiolalaboratory.github.io/sccomp/articles/splines.html) walks through continuous smooths, basis-size choices, prediction and visualisation, the underlying basis decomposition, hierarchical factor smooths, and models containing multiple smooth and random-effect terms. + +We look forward to seeing how the community uses these capabilities to study pseudotime, ageing, treatment response, dose-response relationships, and other continuous biological processes. + +
+Session information + +:::{.smaller} +```{r} +#| label: session-info +#| echo: false + +sessionInfo() +``` +::: + +