4. Case Study: An Inflamed Fibroblast State in Spatial Colitis Data
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).
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 97Step 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
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.3307692
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.
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.
R Session Info
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 26.04.1 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[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: Etc/UTC
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 purrr_1.2.2 xfun_0.61
[25] bluster_1.23.1 cachem_1.1.0 beachmat_2.29.2
[28] jsonlite_2.0.0 DelayedArray_0.39.6 BiocParallel_1.47.0
[31] irlba_2.3.7 parallel_4.6.1 cluster_2.1.8.3
[34] R6_2.6.1 bslib_0.12.0 RColorBrewer_1.1-3
[37] limma_3.99.0 GGally_2.4.0 jquerylib_0.1.4
[40] Rcpp_1.1.2 iterators_1.0.14 knitr_1.52
[43] splines_4.6.1 Matrix_1.7-6 igraph_2.3.3
[46] tidyselect_1.2.1 abind_1.4-8 yaml_2.3.12
[49] viridis_0.6.5 doParallel_1.0.17 codetools_0.2-20
[52] lattice_0.23-1 tibble_3.3.1 withr_3.0.3
[55] S7_0.2.2 evaluate_1.0.5 survival_3.8-12
[58] ggstats_0.14.0 fitdistrplus_1.2-6 circlize_0.4.18
[61] pillar_1.11.1 BiocManager_1.30.27 checkmate_2.3.4
[64] foreach_1.5.2 scales_1.4.0 RhpcBLASctl_0.23-42
[67] glue_1.8.1 metapod_1.21.0 maketools_1.3.2
[70] tools_4.6.1 BiocNeighbors_2.7.3 sys_3.4.3
[73] ScaledMatrix_1.21.0 locfit_1.5-9.12 buildtools_1.0.0
[76] scran_1.41.1 grid_4.6.1 tidyr_1.3.2
[79] colorspace_2.1-3 edgeR_4.99.6 beeswarm_0.4.0
[82] BiocSingular_1.29.1 vipor_0.4.7 cli_3.6.6
[85] rsvd_1.0.5 textshaping_1.0.5 S4Arrays_1.13.1
[88] viridisLite_0.4.3 ComplexHeatmap_2.29.0 dplyr_1.2.1
[91] gtable_0.3.6 isotree_0.6.1-5 sass_0.4.10
[94] digest_0.6.39 SparseArray_1.13.2 ggrepel_0.9.8
[97] dqrng_0.4.1 rjson_0.2.23 farver_2.1.2
[100] htmltools_0.5.9 lifecycle_1.0.5 GlobalOptions_0.1.4
[103] statmod_1.5.2 MASS_7.3-66