--- title: "4. Case Study: An Inflamed Fibroblast State in Spatial Colitis Data" author: - name: Anthony Christidis affiliation: - &core_affiliation Core for Computational Biomedicine, Harvard Medical School - name: Andrew Ghazi affiliation: - *core_affiliation - name: Smriti Chawla affiliation: - *core_affiliation - name: Nitesh Turaga affiliation: - *core_affiliation - name: Ludwig Geistlinger affiliation: - *core_affiliation - name: Robert Gentleman affiliation: - *core_affiliation package: scDiagnostics output: BiocStyle::html_document: toc: true toc_float: true vignette: > %\VignetteIndexEntry{4. Case Study: An Inflamed Fibroblast State in Spatial Colitis Data} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: markdown: wrap: 72 --- ```{r setup, include = FALSE} knitr::knit_hooks$set(pngquant = knitr::hook_pngquant) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", dev = "ragg_png", dpi = 72, fig.retina = 2, fig.align = "center", out.width = "100%", pngquant = "--speed=1 --quality=1-5" ) ``` # Purpose [Vignettes 2](https://ccb-hms.github.io/scDiagnostics/articles/ZeiselBenchmarking.html) and [3](https://ccb-hms.github.io/scDiagnostics/articles/COVIDCaseStudy.html) worked with dissociated scRNA-seq data. `scDiagnostics` operates on `r BiocStyle::Biocpkg("SingleCellExperiment")` objects and does not assume any particular assay - the same diagnostics apply directly to a `r BiocStyle::Biocpkg("SpatialExperiment")` object from imaging-based spatial transcriptomics, without modification, and the same holds for `r BiocStyle::Biocpkg("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). ```{r, message=FALSE} library(scDiagnostics) library(SingleCellExperiment) library(SpatialExperiment) set.seed(100) ``` # 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. ```{r, message=FALSE} data("merfish_reference_data") data("merfish_query_data") class(merfish_reference_data) table(merfish_reference_data$cell_type_merged) table(merfish_query_data$tier2_merged[merfish_query_data$cell_type_merged == "Fibroblast"]) ``` # Step 1: project ```{r, fig.height=6, fig.width=10} 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) ``` # Step 2: detect ```{r, fig.height=5, fig.width=10} 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) ``` ```{r, fig.height=5, fig.width=10} plot(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: ```{r, message=FALSE} fibro_query <- merfish_query_data[, merfish_query_data$cell_type_merged == "Fibroblast"] tapply(is_anomalous, fibro_query$tier2_merged, mean) ``` In 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: ```{r, message=FALSE} 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), ] ``` `plot()` 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: ```{r, fig.height=6, fig.width=8} 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): ```{r, fig.height=5, fig.width=8} 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 SessionInfo, echo=FALSE, message=FALSE, warning=FALSE, comment=NA} options(width = 80) sessionInfo() ```