A few more chnages to day1-1c and started to work on day1-2

Author: julien-roux

assets/images/VisiumHD_nuclei_segmentation_binning_figure1.png

@@ -23,6 +23,7 @@ By the end of this exercise, you will be able to:
 #| warning: false
 #| output: false
 
+library(scuttle)
 library(sf)
 library(arrow)
 library(dplyr)
@@ -33,7 +34,8 @@ library(SpatialExperiment)
 library(SpatialFeatureExperiment)
 library(Voyager)
 library(VisiumIO)
-
+library(ggspavis)
+library(HDF5Array)
 ```
 
 Finally, we will work with the **Visium HD segmented output**. This data corresponds to the same sample and region of interest as the binned data from `Exercise 1A`, but instead of being binned, it contains the output of cell segmentation. Since Space Ranger version 4.0, the output of cell segmentation is available as an optional output of the pipeline, and it contains the coordinates of the segmented cells and their boundaries.
@@ -43,7 +45,7 @@ Finally, we will work with the **Visium HD segmented output**. This data corresp
 ```{r day1-1c-visium_HD_segmented-3}
 options(timeout = 2000)
 
-if (!dir.exists("data/Visium_HD_Human_Colon_Cancer_segmented_outputs/")) {
+if (!dir.exists("data/segmented_outputs/")) {
   dir.create("data", showWarnings = FALSE)
   download.file(
     url = "https://cf.10xgenomics.com/samples/spatial-exp/4.0.1/Visium_HD_Human_Colon_Cancer/Visium_HD_Human_Colon_Cancer_segmented_outputs.tar.gz",
@@ -64,14 +66,21 @@ if (!dir.exists("data/Visium_HD_Human_Colon_Cancer_segmented_outputs/")) {
 
 ::: callout-important
 ## Exercise 1
-Which class do you think should be used to store the Visium HD segmented data? `SpatialExperiment` as we used in Exercise 1A or `SpatialFeatureExperiment` as we used in Exercise 1B?
+- Which class do you think should be used to store the Visium HD segmented data? `SpatialExperiment` as we used in Exercise 1A or `SpatialFeatureExperiment` as we used in Exercise 1B?
+- What does the data for each segmented cell represent? What is the binning strategy? What are the potential limitations?
+
 ::: 
 
 ::: {.callout-tip collapse="true"}
 ## Answer
-Since segmented data contains cell boundaries, the coordinates of the segmented cells are not limited to the spot coordinates as in the binned data. Therefore, we will use `SpatialFeatureExperiment` to work with the segmented data, as it allows us to store and visualize spatial features such as cell boundaries, as we did for the Xenium data in Exercise 1B. 
+- Since segmented data contains cell boundaries, the coordinates of the segmented cells are not limited to the spot coordinates as in the binned data. Therefore, we will use `SpatialFeatureExperiment` to work with the segmented data, as it allows us to store and visualize spatial features such as cell boundaries, as we did for the Xenium data in Exercise 1B. 
+`SpatialExperiment` is more suitable for binned data where the coordinates correspond to the spot centroids on a regular grid and there are no additional spatial features to store.
+
+- The Visium HD still captures the data on 2x2 µm bins. These are partitioned into larger bins based on the cell they correspond to (see Figure below). The partitioning of aggregated bins mimics single-cell data since the UMI counts are reported per segmented cell basis.
+
+![Binning approach in the Visium HD segmented output data (segmented nuclei here)](../assets/images/VisiumHD_nuclei_segmentation_binning_figure1.png)
 
-`SpatialExperiment` is more suitable for binned data where the coordinates correspond to the spot centroids and there are no additional spatial features to store.
+One limitation is that a bin can overlap than one cell, so the assignment of the transcripts captured would be ambiguous. Of course this approach relies also on the quality of the segmentation mask. 
 ::: 
 
 
@@ -123,7 +132,7 @@ spe <- loadHDF5SummarizedExperiment(dir="results/day1/", prefix="01.1_spe_")
 ```
 
 ::: callout-important
-## Exercise 1
+## Exercise 2
 
 - 1. Compare the number of features and spots in the segmented and binned objects. What do you observe?
 - 2. Compare the `colData` of the two objects. What are the differences in the metadata columns? What do they represent?
@@ -196,15 +205,65 @@ plotSpatialFeature(sfe_seg,
                    features="PIGR", 
                    exprs_values = "counts") 
 ```  
+:::
+
+::: callout-important
+## Exercise 3
+Plot the cell segmentation mask of very small region of interest on the slide, for example:
 
-Save the object
 ```{r day1-1c-visium_HD_segmented-10}
+roi_zoomed <- c(
+  xmin = 54500,
+  xmax = 55000,
+  ymin = 12500,
+  ymax = 13000
+)
+```
+What do you observe?
+::: 
+
+::: {.callout-tip collapse="true"}
+## Answer
+
+```{r day1-1c-visium_HD_segmented-11}
+sfe_zoomed <- sfe_seg[, spatialCoords(sfe_seg)[, 1] >= roi_zoomed[["xmin"]] &
+                        spatialCoords(sfe_seg)[, 1] <= roi_zoomed[["xmax"]] &
+                        spatialCoords(sfe_seg)[, 2] >= roi_zoomed[["ymin"]] &
+                        spatialCoords(sfe_seg)[, 2] <= roi_zoomed[["ymax"]]]
+
+plot(st_geometry(colGeometries(sfe_zoomed)[["cellSeg"]]), 
+     col=rep(colors(), 2), 
+     lty=0)
+rm(sfe_zoomed)
+```
+
+When zooming in, we can clearly see that 2 um bins make up each cell in the dataset.
+:::
+
+
+Save the object
+```{r day1-1c-visium_HD_segmented-12}
 saveHDF5SummarizedExperiment(sfe_seg,
   dir = "results/day1", prefix = "01.1_sfe_visium_", replace = TRUE,
   chunkdim = NULL, level = NULL, as.sparse = NA,
   verbose = NA
 )
 ```
 
-For the rest of the exercises, we will use this Visium HD segmented object, but feel free to replace by the binned object or the Xenium object to compare the results.
+<!--
+TO DO: Even using this, some out-of-memory components break when re-importing because the original file is unavailable or is tied to some absolute path. E.g., 
+
+Error in (function (cond)  : 
+  error in evaluating the argument 'x' in selecting a method for function 'colData': external pointer is not valid
+
+Another (better) way to save the object would be to use the alabaster.sfe package, but I failed to install it on my Mac...
+
+```{r}
+fsave <- file.path("results/day1", "01.1_sfe_visium")
+alabaster.sfe::saveObject(sfe_seg, fsave)
+dir_tree(fsave)
+```
+-->
+
+For the rest of the exercises, we will use this Visium HD segmented object, but feel free to replace by the binned object (`Exercise 1A`) or the Xenium object (`Exercise 1B`) to compare the results.
 

day1/day1-1c_visium_HD_segmented.qmd

@@ -9,16 +9,17 @@ execute:
 ---
 
 
-## Quality control at bin-level
+## Quality control at cell-level
 
-In this second exercise, we will focus on the critical step of quality control (QC) for sequencing-based spatial transcriptomics data. The Visium HD object contains 16 um bins, so each observation is a spatial bin rather than an individual cell. QC helps flag low-quality bins or technical artifacts before downstream analysis.
+In this second exercise, we will focus on the critical step of quality control (QC) for sequencing-based spatial transcriptomics data. We will focus on the segmented Visium HD dataset used in `Exercise 1C`, so each observation is an individual cell with the selected region of interest. QC helps flag low-quality cells or technical artifacts before downstream analysis.
 
+*NOTE:* A similar QC can be done per bin or aggregated bin if a binned analysis of Visium HD is performed.
 
 ## Learning objectives
 
 By the end of this exercise, you will be able to:
 
-- Calculate per-bin QC metrics.
+- Calculate per-cell QC metrics.
 - Identify local outliers based on various QC metrics.
 - Detect spatial artifacts using `SpotSweeper`.
 - Visualize QC metrics and detected artifacts.
@@ -31,41 +32,45 @@ By the end of this exercise, you will be able to:
 #| output: false
 
 library(SpatialExperiment)
+library(SpatialFeatureExperiment)
+library(Voyager)
+library(HDF5Array)
+
 library(SpotSweeper)
 library(scuttle)
+library(scrapper)
 library(scater)
-library(ggside)
-library(ggplot2)
 library(escheR) 
-library(HDF5Array)
+
+library(ggplot2)
+library(ggside)
 library(ggspavis)
 library(patchwork)
 ```
 
 ## Calculate QC metrics
 
-We will start by calculating QC metrics that are also commonly used in scRNA-seq data analysis: total UMI counts, number of detected genes, and percentage of UMIs originating from mitochondrial genes. Here these metrics are calculated per Visium HD bin and used to identify low-quality bins.
+We will start by calculating QC metrics that are also commonly used in scRNA-seq data analysis: total UMI counts, number of detected genes, and percentage of UMIs originating from genes of the mitochondrial genome. Here these metrics are calculated per segmented cell and used to identify low-quality cells.
 
-We will start by loading our `SpatialExperiment` object, which was saved in the previous exercise, and prepare it for quality control analysis.  
+We will start by loading our `SpatialFeatureExperiment` object, which was saved in the previous exercise, and prepare it for quality control analysis.  
 
 
 ```{r}
-# Load the SpatialExperiment object
-sfe <- loadHDF5SummarizedExperiment(dir="results/day1/", prefix="01.1_sfe_seg_")
+# Load the SpatialExperiment object if not in the R session already
+sfe <- loadHDF5SummarizedExperiment(dir="results/day1/", prefix="01.1_sfe_visium_")
 
 # Identify mitochondrial genes (genes with symbol starting with "MT-")
 is.mito <- rownames(sfe)[grepl("^MT-", rownames(sfe))]
 ```
 
 <!--- NOTE: HDF5-backed saving/loading via saveHDF5SummarizedExperiment() is reliable for assays, but it can leave sf geometries with invalid external pointers after reload (GEOS/s2), which then crash functions like localOutliers() with “external pointer is not valid”. 
-Need a solution or don't delete previous object --->
+Need a solution or don't delete previous object 
+--->
 
-Now, we will calculate per-bin QC metrics using the `scuttle` package. The function `addPerCellQCMetrics` adds QC metrics, including mitochondrial gene percentage, to the `colData` of the `SpatialExperiment` object. The function name says "cell" because it comes from single-cell workflows, but here it operates on spatial bins.
-
-Run the following code to compute these metrics and inspect the results. 
+Now, we will calculate per-cell QC metrics using the `scrapper` package. The function `quickRnaQc.se` adds QC metrics, including mitochondrial gene percentage, to the `colData` of the `SpatialFeatureExperiment` object. Run the following code to compute these metrics and inspect the results:
 
 ```{r}
-sfe <- addPerCellQCMetrics(sfe, subsets = list(Mito = is.mito))
+sfe <- quickRnaQc.se(sfe, subsets = list(Mito = is.mito))
 ```
 
 ::: callout-important
@@ -82,15 +87,14 @@ Check which metadata has been added to `colData`. Do you recognize the different
 colData(sfe)
 ```
 
-- The `sum` column contains the total number of unique molecular identifiers (UMIs) for each bin
-- The `detected` column contains the number of unique genes detected per bin
-- The `subsets_Mito_percent` column contains the percentage of transcripts mapping to mitochondrial genes per bin
+- The `sum` column contains the total number of unique molecular identifiers (UMIs) for each cell
+- The `detected` column contains the number of unique genes detected per cell
+- The `subset.proportion.Mito` column contains the percentage of transcripts mapping to mitochondrial genes per cell
 :::
 
-
 ### Global outlier detection
 
-In order to detect low-quality bins and filter them out, QC methods adapted from scRNA-seq analysis can be applied. A simple option is to apply a fixed upper or lower threshold to a metric across all bins and remove the bins that do not pass the filtering criteria.
+In order to detect low-quality cells and filter them out, QC methods adapted from scRNA-seq analysis can be applied. A simple option is to apply a fixed upper or lower threshold to a metric across all cells and remove the cells that do not pass the filtering criteria.
 
 To set up these thresholds we can look for information from 10x Genomics and the community, for example https://github.com/10XGenomics/HumanColonCancer_VisiumHD/issues/28. Based on this discussion the authors of the OSTA book wrote (see "Remove bins overlaying empty tissue" section in https://lmweber.org/OSTA/pages/seq-workflow-visium-hd-bin.html#dependencies):
 
@@ -127,15 +131,15 @@ p2 <- ggcells(sfe, aes(x=detected, y=subsets_Mito_percent)) +
 p1 + p2
 ```
 
-And further check how many bins would be removed with the given thresholds:
+And further check how many cells would be removed with the given thresholds:
 
 ```{r}
-# Flag bins based on fixed thresholds
+# Flag cells based on fixed thresholds
 sfe$qc_lib_size <- sfe$sum < 400
 sfe$qc_detected <- sfe$detected < 300
 sfe$qc_mito_prop <- sfe$subsets_Mito_percent > 20
 
-# Tabulate number of bins flagged by each
+# Tabulate number of cells flagged by each
 qc <- grep("^qc", names(colData(sfe)))
 sapply(colData(sfe)[qc], table)
 ```
@@ -144,7 +148,7 @@ sapply(colData(sfe)[qc], table)
 ## Exercise 2
 What do you think about this initial filtering?
 
-Have a look at the mitochondrial percentage values on the tissue section using `plotVisium()`. Which bins are more likely to be excluded?
+Have a look at the mitochondrial percentage values on the tissue section using `plotVisium()`. Which cells are more likely to be excluded?
 :::
 
 :::{.callout-tip collapse="true"}
@@ -165,7 +169,7 @@ plot_layout(nrow=1) &
 Very importantly, we want to have a look into the spatial pattern of using a fixed threshold for QC metrics :
 
 ```{r}
-# check spatial pattern of discarded bins
+# check spatial pattern of discarded cells
 p1 <- plotObsQC(sfe, 
                 plot_type="spot", annotate="qc_lib_size") + 
   ggtitle("Library size (< 400 UMI)")
@@ -185,48 +189,48 @@ gridExtra::grid.arrange(p1, p2, p3, ncol=2)
 
 ### Local Outlier Detection
 
-To avoid discarding whole layers of tissue, another approach to identify low-quality bins is to use local outlier detection methods that take into account the spatial context of each bin.
+To avoid discarding whole layers of tissue, another approach to identify low-quality cells is to use local outlier detection methods that take into account the spatial context of each cell.
 
 `SpotSweeper` package provides functions to identify common QC metrics (such as the ones we have seen in the previous section: library size, number of detected genes, and mitochondrial percentage) and detect local outliers based on these metrics.
 
-We will use `localOutliers` from `SpotSweeper`, which helps identify bins that deviate significantly from their local neighborhood.
+We will use `localOutliers` from `SpotSweeper`, which helps identify cells that deviate significantly from their local neighborhood.
 
 
 ```{r}
 # Identify local outliers based on library size ("sum" of counts).
-# Bins with unusually low library size compared to their neighbors will be flagged.
+# cells with unusually low library size compared to their neighbors will be flagged.
 sfe <- localOutliers(sfe,
   metric = "sum",
   direction = "lower",
   log = TRUE
 )
 
 # Identify local outliers based on the number of unique genes detected ("detected").
-# Bins with an unusually low number of detected genes will be flagged.
+# cells with an unusually low number of detected genes will be flagged.
 sfe <- localOutliers(sfe, 
   metric = "detected",
   direction = "lower",
   log = TRUE
 )
 
 # Identify local outliers based on the mitochondrial gene percentage.
-# Bins with an unusually high mitochondrial percentage will be flagged.
+# cells with an unusually high mitochondrial percentage will be flagged.
 sfe <- localOutliers(sfe,
   metric = "subsets_Mito_percent",
   direction = "higher",
   log = FALSE
 )
 
 # Combine all individual outlier flags into a single "local_outliers" column.
-# A bin is considered a local outlier if it is flagged by any of the above metrics.
+# A cell is considered a local outlier if it is flagged by any of the above metrics.
 sfe$local_outliers <- as.logical(sfe$sum_outliers) |
   as.logical(sfe$detected_outliers) |
   as.logical(sfe$subsets_Mito_percent_outliers)
 ```
 
 ::: callout-important
 ## Exercise 3
-How many bins were identified as local outliers based on the combined criteria?
+How many cells were identified as local outliers based on the combined criteria?
 
 As we have seen before, visualizing the QC metrics is crucial for understanding the quality of your spatial transcriptomics data and for deciding on appropriate filtering strategies, so look at their spatial localization using the `plotObsQC` function.
 
@@ -237,7 +241,7 @@ What do you think?
 :::{.callout-tip collapse="true"}
 ### Answer
 ```{r}
-# Count the number of bins flagged as local outliers.
+# Count the number of cells flagged as local outliers.
 table(sfe$local_outliers)
 
 plotObsQC(sfe, 
@@ -283,7 +287,7 @@ plotQCmetrics(sfe,
   ggtitle("Local outliers for number of UMIs")
 ```
 
-Very few bins seem to be flagged as outliers!
+Very few cells seem to be flagged as outliers!
 :::
 
 
@@ -323,16 +327,16 @@ sfe$discard <- sfe$global_outliers | sfe$local_outliers
 
 ::: callout-important
 ## Exercise 5
-Visualize the bins that will be discarded based on the combined criteria. How many bins will be removed?
+Visualize the cells that will be discarded based on the combined criteria. How many cells will be removed?
 :::
 
 :::{.callout-tip collapse="true"}
 ### Answer
 ```{r}
-# check how many bins will be discarded
+# check how many cells will be discarded
 table(sfe$discard)
 
-# have a final look at the bins that we will discard
+# have a final look at the cells that we will discard
 plotSpatialFeature(sfe, colGeometryName = "cellseg", features = "discard")
 ```
 :::
@@ -388,10 +392,10 @@ plotQCmetrics(sfe,
 :::
 
 ## Filtering 
-We decide to filter out the low-quality bins from the `SpatialExperiment` object based on the combined criteria of global thresholds and local outliers. Additionally, we will remove any features (genes) that have zero counts across all remaining bins.
+We decide to filter out the low-quality cells from the `SpatialExperiment` object based on the combined criteria of global thresholds and local outliers. Additionally, we will remove any features (genes) that have zero counts across all remaining cells.
 
 ```{r}
-# remove the bins considered of low quality based on both global and local metrics
+# remove the cells considered of low quality based on both global and local metrics
 sfe <- sfe[, !sfe$discard]
 
 # remove features with all 0 counts
@@ -421,7 +425,7 @@ gc()
 **Key Takeaways:**
 
 -   `scuttle::addPerCellQCMetrics` provides essential per-bin QC metrics.
--   `SpotSweeper::localOutliers` helps identify individual problematic bins.
+-   `SpotSweeper::localOutliers` helps identify individual problematic cells.
 -   `SpotSweeper::findArtifacts` is powerful for detecting spatially coherent regions of poor quality.
 -   Visualizing QC results is critical for data interpretation and downstream analysis decisions.
 :::