4. Case Study: An Inflamed Fibroblast State in Spatial Colitis Data
23 September 2026
Package
scDiagnostics 1.7.15
1 Purpose
Vignettes 2
and 3
worked with dissociated scRNA-seq data. scDiagnostics operates on
SingleCellExperiment objects and does not
assume any particular assay - the same diagnostics apply directly to a
SpatialExperiment object from imaging-based
spatial transcriptomics, without modification, and the same holds for
SpatialFeatureExperiment objects (which extend
SpatialExperiment with explicit cell/tissue geometries), since both
inherit the same core SingleCellExperiment interface that
scDiagnostics relies on. This vignette repeats the project-detect-
characterize workflow on MERFISH spatial data from a mouse model of
DSS-induced colitis (Cadinu et al. 2024): healthy colon tissue (Day 0,
reference) versus tissue at peak inflammation (Day 9, query).
library(scDiagnostics)
library(SingleCellExperiment)
library(SpatialExperiment)
set.seed(100)2 The data
merfish_reference_data and merfish_query_data are downsampled
subsets on the same 943-gene targeted MERFISH panel, restricted to 5
shared cell types; see ?merfish_reference_data for details. At Day 9,
some fibroblasts further split into an inflammation-associated
“Inflamed Fibroblast” state (tier2_merged) that does not exist at Day
0; a merged cell_type_merged column collapses this back into a shared
“Fibroblast” label so the two timepoints can be compared directly.
data("merfish_reference_data")
data("merfish_query_data")
class(merfish_reference_data)
#> [1] "SpatialExperiment"
#> attr(,"package")
#> [1] "SpatialExperiment"
table(merfish_reference_data$cell_type_merged)
#>
#> Endothelial Epithelial Fibroblast Other Immune Smooth Muscle
#> 220 220 220 220 220
table(merfish_query_data$tier2_merged[merfish_query_data$cell_type_merged == "Fibroblast"])
#>
#> Fibroblast Inflamed Fibroblast
#> 163 973 Step 1: project
plotCellTypePCA(
query_data = merfish_query_data,
reference_data = merfish_reference_data,
cell_types = unique(merfish_reference_data$cell_type_merged),
query_cell_type_col = "cell_type_merged",
ref_cell_type_col = "cell_type_merged",
pc_subset = 1:3)
#> Picking joint bandwidth of 0.279
#> Picking joint bandwidth of 0.222
#> Picking joint bandwidth of 0.251
4 Step 2: detect
anomaly_output <- detectAnomaly(
reference_data = merfish_reference_data,
query_data = merfish_query_data,
ref_cell_type_col = "cell_type_merged",
query_cell_type_col = "cell_type_merged",
cell_types = "Fibroblast",
pc_subset = 1:5,
n_tree = 500)
is_anomalous <- anomaly_output[["Fibroblast"]]$query_anomaly
mean(is_anomalous)
#> [1] 0.3307692plot(anomaly_output, cell_type = "Fibroblast", pc_subset = 1:3, data_type = "query")
Because merfish_query_data retains the original fine-grained
tier2_merged label, we can check how the flagged cells relate to the
ground-truth “Inflamed Fibroblast” state:
fibro_query <- merfish_query_data[, merfish_query_data$cell_type_merged == "Fibroblast"]
tapply(is_anomalous, fibro_query$tier2_merged, mean)
#> Fibroblast Inflamed Fibroblast
#> 0.1779141 0.5876289In this downsampled dataset, cells with the ground-truth “Inflamed Fibroblast” label are flagged as anomalous roughly twice as often as plain “Fibroblast” cells - a real enrichment, though far from a clean separation, consistent with inflammation being a graded rather than binary state at the single-cell level.
5 Step 3: characterize
To characterize what distinguishes the flagged cells, calculateGeneShifts()
can run its own internal anomaly detection (detect_anomalies = TRUE,
anomaly_comparison = TRUE) on the full 943-gene panel - the same
panel detectAnomaly() used above - while restricting the actual
statistical comparison to a small extracellular matrix (ECM) gene panel
via genes_to_analyze. This keeps detection and characterization
consistent without needing to manually subset cells first:
ecm_signature <- c("Col1a2", "Timp2", "Col6a1", "Sparc", "Dpt")
gene_shifts <- calculateGeneShifts(
query_data = merfish_query_data,
reference_data = merfish_reference_data,
query_cell_type_col = "cell_type_merged",
ref_cell_type_col = "cell_type_merged",
cell_types = "Fibroblast",
pc_subset = 1:5,
genes_to_analyze = ecm_signature,
detect_anomalies = TRUE,
anomaly_comparison = TRUE)
gene_shifts$PC1[order(gene_shifts$PC1$p_adjusted), ]
#> gene loading cell_type p_value mean_query mean_reference p_adjusted
#> 1 Dpt NA Fibroblast 5.122401e-26 0.1317246 1.991780 2.561201e-25
#> 2 Sparc NA Fibroblast 5.230999e-23 0.2882483 1.868690 1.307750e-22
#> 3 Col1a2 NA Fibroblast 1.030939e-22 0.4023466 2.106959 1.718231e-22
#> 4 Timp2 NA Fibroblast 9.084576e-20 0.1081913 1.289465 1.135572e-19
#> 5 Col6a1 NA Fibroblast 2.190372e-18 0.1925289 1.396044 2.190372e-18
#> significant
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 TRUE
#> 5 TRUEplot() shows this as a heatmap of per-gene z-scores, one column per
cell (reference and query, annotated by anomaly status), rather than
collapsing each group to a single averaged column - which matters here,
since it reveals more than the group averages alone would:
plot(gene_shifts, cell_type = "Fibroblast", pc_subset = 1:5,
plot_type = "heatmap", plot_by = "p_adjusted", n_genes = 5,
show_anomalies = TRUE)
All five ECM genes are significantly lower in the anomalous query
fibroblasts than in the reference, but the heatmap shows this isn’t a
single uniform effect. Col1a2 and Sparc show a visible extra drop
concentrated specifically in the anomalous (rightmost, red-annotated)
cells, beyond what’s already true of the query more broadly. Timp2,
Col6a1, and Dpt, on the other hand, are already substantially
reduced across essentially all query fibroblasts - anomalous or not -
so for those three genes, detectAnomaly() isn’t isolating a
distinctly-shifted subgroup so much as reflecting a shift already
present dataset-wide. This split held up across several PC-subset
choices we checked, so it looks like a real property of this dataset
rather than a detection parameter to tune away: some genes distinguish
the specific cells flagged as anomalous, and others distinguish the
query condition as a whole - both are useful, but they’re different
claims.
The barplot below focuses on Col1a2 and Sparc specifically, since
they’re the clearest example of the anomaly-specific pattern - a
“Query Non-Anomaly vs Ref” bar close to zero alongside a much larger
“Query Anomaly vs Ref” bar, showing the shift really is concentrated in
the flagged cells for these two genes (unlike Timp2/Col6a1/Dpt
above, where the non-anomaly bar would already be nearly as large as
the anomaly bar):
gene_shifts_focused <- calculateGeneShifts(
query_data = merfish_query_data,
reference_data = merfish_reference_data,
query_cell_type_col = "cell_type_merged",
ref_cell_type_col = "cell_type_merged",
cell_types = "Fibroblast",
pc_subset = 1:5,
genes_to_analyze = c("Col1a2", "Sparc"),
detect_anomalies = TRUE,
anomaly_comparison = TRUE)
plot(gene_shifts_focused, cell_type = "Fibroblast", pc_subset = 1:5,
plot_type = "barplot", plot_by = "p_adjusted", n_genes = 2,
show_anomalies = TRUE, pseudo_bulk = TRUE)
As with the COVID-19 case study, this is a specific finding about this cell population in this dataset; the broader claim that this workflow generalizes across data modalities is best supported by comparing this result to the scRNA-seq case study in vignette 3, not by either result alone.
6 R Session Info
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.5 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB LC_COLLATE=C
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] SpatialExperiment_1.23.0 SingleR_2.15.4
[3] scater_1.41.2 ggplot2_4.0.3
[5] scuttle_1.23.2 splatter_1.37.1
[7] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
[9] Biobase_2.73.2 GenomicRanges_1.65.4
[11] Seqinfo_1.3.2 IRanges_2.47.5
[13] S4Vectors_0.51.10 BiocGenerics_0.59.12
[15] generics_0.1.4 MatrixGenerics_1.25.0
[17] matrixStats_1.5.0 scDiagnostics_1.7.15
[19] BiocStyle_2.41.0
loaded via a namespace (and not attached):
[1] gridExtra_2.3.1 rlang_1.3.0 magrittr_2.0.5
[4] clue_0.3-68 GetoptLong_1.1.1 otel_0.2.0
[7] ggridges_0.5.7 compiler_4.6.1 png_0.1-9
[10] systemfonts_1.3.2 vctrs_0.7.3 shape_1.4.6.1
[13] crayon_1.5.3 pkgconfig_2.0.3 fastmap_1.2.0
[16] backports_1.5.1 magick_2.9.1 XVector_0.53.0
[19] labeling_0.4.3 rmarkdown_2.32 ggbeeswarm_0.7.3
[22] ragg_1.5.2 tinytex_0.61 purrr_1.2.2
[25] xfun_0.61 bluster_1.23.1 cachem_1.1.0
[28] beachmat_2.29.3 jsonlite_2.0.0 DelayedArray_0.39.7
[31] BiocParallel_1.47.0 irlba_2.3.7 parallel_4.6.1
[34] cluster_2.1.8.3 R6_2.6.1 bslib_0.12.0
[37] RColorBrewer_1.1-3 limma_3.99.0 GGally_2.4.0
[40] jquerylib_0.1.4 iterators_1.0.14 Rcpp_1.1.2
[43] bookdown_0.48 knitr_1.52 splines_4.6.1
[46] Matrix_1.7-6 igraph_2.3.3 tidyselect_1.2.1
[49] dichromat_2.0-1 abind_1.4-8 yaml_2.3.12
[52] viridis_0.6.5 doParallel_1.0.17 codetools_0.2-20
[55] lattice_0.23-1 tibble_3.3.1 withr_3.0.3
[58] S7_0.2.2 evaluate_1.0.5 survival_3.8-12
[61] ggstats_0.14.0 fitdistrplus_1.2-6 circlize_0.4.18
[64] pillar_1.11.1 BiocManager_1.30.27 checkmate_2.3.4
[67] foreach_1.5.2 scales_1.4.0 RhpcBLASctl_0.23-42
[70] glue_1.8.1 metapod_1.21.0 tools_4.6.1
[73] BiocNeighbors_2.7.3 ScaledMatrix_1.21.0 locfit_1.5-9.12
[76] scran_1.41.1 Cairo_1.7-0 grid_4.6.1
[79] tidyr_1.3.2 colorspace_2.1-3 edgeR_4.99.6
[82] beeswarm_0.4.0 BiocSingular_1.29.1 vipor_0.4.7
[85] cli_3.6.6 rsvd_1.0.5 textshaping_1.0.5
[88] S4Arrays_1.13.1 viridisLite_0.4.3 ComplexHeatmap_2.29.0
[91] dplyr_1.2.1 gtable_0.3.6 isotree_0.6.1-5
[94] sass_0.4.10 digest_0.6.39 SparseArray_1.13.3
[97] ggrepel_0.9.8 dqrng_0.4.1 rjson_0.2.23
[100] farver_2.1.2 htmltools_0.5.9 lifecycle_1.0.5
[103] GlobalOptions_0.1.4 statmod_1.5.2 MASS_7.3-66