--- title: "Comparing proximity effects between conditions" author: "Elijah Willie" date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true toc_float: true vignette: > %\VignetteIndexEntry{Comparing proximity effects between conditions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, fig.width = 6, fig.height = 4, dpi = 90) ``` # Introduction The `r Biocpkg("PACE")` manual fits one section. This vignette asks whether a proximity effect differs between two groups of patients. That needs a cohort, and adds one block to the model. Beside the spatial cell state block, the proximity effect common to both arms, PACE estimates a **responder spatial state** block: the part that differs between them. # The melanoma subset The package ships a subset of a CosMx melanoma cohort of patients profiled after immune checkpoint blockade, contrasting progressive disease (PD) against non-progressive disease (non-PD, pooling stable, partial and complete responders). ```{r load, message = FALSE} library(PACE) library(SpatialExperiment) library(ggplot2) library(dplyr) spe <- readRDS(system.file("extdata", "mel_cosmx_subset.rds", package = "PACE")) spe table(spe$cellType) table(unique(data.frame(image = spe$image, arm = spe$Responder))$arm) ``` One section per patient, so the contrast is between patients. Every cell and patient of the published cohort is here; only the panel is reduced, to the 180 most detected of 927 genes. Derivation: `inst/scripts/make-mel-cosmx-subset.R`. > **Provenance.** These data are a subset of the deposit accompanying Dong, Su, > Kluger, Fan and Kluger, *SIMVI disentangles intrinsic and spatial-induced > cellular states in spatial omics data* > ([doi:10.5281/zenodo.14708000](https://doi.org/10.5281/zenodo.14708000)), > used here under CC-BY-4.0. # Fitting with a condition Three arguments beyond the single-section call: - `condition_col` names the contrast and switches on the responder block. - `kernel_per_image = TRUE` keeps neighbourhoods within a section, so no cell neighbours another patient's. - `image_re = "intercept"` absorbs patient-level shifts in overall expression. ```{r fit} fit <- paceFit(spe, celltype_col = "cellType", condition_col = "Responder", image_col = "image", kernel_per_image = TRUE, image_re = "intercept", dispersion = "nb1", verbose = FALSE) fit ``` The interaction term is named from the non-reference level of `condition_col`: ```{r respterm} paceParams(fit)$resp_term ``` so `ResponderPD:Tumour` is the extra change in expression per unit of tumour proximity in the PD arm, on top of the `Tumour` effect shared by both arms. # The responder block in the decomposition `varianceDecomposition()` now carries a `Responder spatial %` column beside the `Spatial %` one: ```{r decomp} varianceDecomposition(fit) |> group_by(focal) |> summarise(`Spatial %` = round(median(`Spatial %`), 4), `Responder spatial %` = round(median(`Responder spatial %`), 4), `Spillover %` = round(median(`Spillover %`), 2), .groups = "drop") ``` The responder block is much the smaller: most of a cell's response to its neighbours is common to both arms. Contamination exceeds both, and is largest in the sparse immune populations. `plotDecomposition()` shows the same thing. On a fit with a condition it adds a "Responder spatial state" block beside the shared spatial one, so the two are directly comparable; on a fit without one that block is simply absent. The right panel zooms in on the small components, which are otherwise invisible next to cell-type identity. ```{r decomp-plot, fig.width = 11, fig.height = 4.5} plotDecomposition(fit) ``` # Which pairs carry the condition difference `pairVariance()` and `plotPairHeatmap()` both take `block = "responder"`, which attributes the responder block across neighbours instead of the shared spatial one: ```{r pairs-responder, fig.width = 6, fig.height = 4.6} pairVariance(fit, block = "responder") |> arrange(desc(val)) |> head(5) plotPairHeatmap(fit, block = "responder") ``` Fibroblast next to endothelium takes the largest share. Compare the shared spatial block, which orders differently: ```{r pairs-spatial, fig.width = 6, fig.height = 4.6} plotPairHeatmap(fit, block = "spatial") ``` # The genes that carry it Pair shares say where to look, `neighbourSlopes()` which genes. The `ResponderPD:` terms differ between arms: ```{r calls} ns <- neighbourSlopes(fit) resp <- subset(ns, grepl("^ResponderPD:", term) & lfsr < 0.05) nrow(resp) resp[order(resp$lfsr), c("gene", "focal", "neighbour", "estimate_shrunk", "lfsr")] |> head(8) ``` `plotDrivers()` gives a pair as a figure: the genes ranked by driver score on the left, their per-gene decomposition on the right. On a fit with a condition the score is built from the `resp_term` interaction, so this ranks genes by how much the proximity effect **differs between arms**, not by the effect they share. Taking tumour cells next to T cells, *CRYAB*, *MZT2A* and *S100B* lead, and nine of the thirteen genes called are MHC class I or interferon: *B2M*, *HLA-B*, *HLA-C*, *HLA-E*, *STAT1*, *MX1*, *IFITM1*, *IFITM3* and *IFI27*. Those nine all have positive interactions, so their response to T cell proximity is steeper in the progressive arm. ```{r drivers, fig.width = 10, fig.height = 4.5} plotDrivers(fit, "Tumour", "T_Cell") ``` Two things follow from the score being the interaction. A difference between arms subtracts anything common to both, and contamination from the neighbour is largely common to both, so this ranking is less exposed to it than a ranking of the shared effect would be. That matters here, because *B2M* and *HLA-B* are expressed by T cells themselves and the right panel shows spillover as their largest block; it is the differencing, not the size of that block, that makes the ranking readable. To rank the *shared* effect instead, re-run `paceShrink()` with the responder term switched off, which puts the shrinkage on the baseline slopes. No new model fit is needed. `topDrivers()` returns the same ranking as a table, with the score and the parts it is built from: ```{r drivers-table} topDrivers(fit)[["Tumour_T_Cell"]]$scores |> select(rank, gene, driver_score, b_clean, spec, focal_mean, lfsr) |> head(8) ``` The table is also filtered at `lfsr < 0.05`, so a sparse figure means few genes cleared that threshold rather than few genes being involved. On this reduced panel `Macrophage <- Tumour` has only one, which is why the pair shown here is a different one. Among the strongest is *SPP1* in macrophages next to tumour cells: ```{r spp1} subset(ns, gene == "SPP1" & focal == "Macrophage" & neighbour == "Tumour") ``` `Tumour` is the response shared by both arms; `ResponderPD:Tumour` is how much it differs in progressive disease. Negative, so macrophage *SPP1* rises less steeply with tumour proximity when disease progressed. This is the paper's melanoma result. The counts behind it, split by arm: ```{r proximity, fig.width = 7, fig.height = 4} plotProximity(fit, spe, "SPP1", "Macrophage", "Tumour", breaks = c(0, 10, 20, 30, 40, Inf), condition = TRUE) ``` `plotResponseCurve()` puts the model over the data: solid lines are binned means and standard errors, dashed lines the PACE slopes read from the fit. ```{r curve, fig.width = 7, fig.height = 5} plotResponseCurve(fit, spe, "SPP1", "Macrophage", "Tumour") ``` The dashed lines diverge; that divergence is the interaction. The x axis is the model's own covariate, so the slopes are drawn on the scale they were estimated on. `plotResponseMap()` gives the spatial view, using the two sections shown in the paper: ```{r map, fig.width = 9, fig.height = 8} plotResponseMap(fit, spe, "SPP1", "Macrophage", "Tumour", images = c(nonPD = "32157_18", PD = "32156_17")) ``` # Session info ```{r sessioninfo} sessionInfo() ```