---

title: "Preprocessing"

output:

  html_document:

    df_print: paged

editor_options: 

  chunk_output_type: console

---

# Preprocessing Data (Following OSTA Analysis Pipeline)

In this markdown, analysis and quality control steps are performed following chapters 8 and 9 of the book 'Orchestrating Spatially-Resolved Transcriptomics Analysis with Bioconductor' <https://lmweber.org/OSTA-book/>. Data is fetched and plotted, then QC metrics are calculated and data is filtered based on QC thresholds for a single sample.

### Overview of Data

Data is sourced from spatialLIBD (<http://spatial.libd.org/spatialLIBD/>), a project containing LIBD human dorsolateral pre-frontal cortex (DLPFC) spatial transcriptomics data from the 10x Genomics Visium Platform (<https://10xgenomics.com/products/spatial-gene-expression>). The data contains 12 lowres samples each with a scale factor of 0.0450045.

### Importing Libraries

```{r Import Libraries, message=FALSE}

library(spatialLIBD)

library(ggspavis)

library(scater)

library(ggplot2)

```

### Helper Functions for QC

The following functions are called by the driver code to perform plotting and QC operations.

```{r Helper Functions for QC, message=FALSE}

#' Plots a SpatialExperiment object by transforming the image into a rasterGrob object and calling ggplot

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @return None

plot_image <- function(spe_sub) {

  # get histology image

  img <- SpatialExperiment::imgRaster(spe_sub)

  ## Transform to a rasterGrob object

  grob <- grid::rasterGrob(img, width = grid::unit(1, "npc"), height = grid::unit(1, "npc"))

  # Create df for plotting

  sample_df <- as.data.frame(colData(spe_sub), optional = TRUE)

  ## Make a plot using geom_spatial

  p <- ggplot2::ggplot(

      sample_df,

      ggplot2::aes(

          x = pxl_col_in_fullres * SpatialExperiment::scaleFactors(spe_sub),

          y = pxl_row_in_fullres * SpatialExperiment::scaleFactors(spe_sub),

      )

  ) +

      geom_spatial(

          data = tibble::tibble(grob = list(grob)),

          ggplot2::aes(grob = grob),

          x = 0.5,

          y = 0.5

      )

  ## Show the plot

  print(p)

  p

}

#' Identifies mitochondrial genes, and adds per spot QC metrics including mitochondrial genes as a feature subset

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @return spe_sub Spatial Experiment object with added QC metrics

add_QC_metrics <- function(spe_sub) {

  # identify mitochondrial genes

  is_mito <- grepl("(^MT-)|(^mt-)", rowData(spe_sub)$gene_name)

  # add per spot QC metrics, including mitochondrial genes as a feature subset

  spe_sub <- addPerCellQC(spe_sub, subsets = list(mito = is_mito))

  spe_sub

}

#' Sets a lower filtering threshold for library size, plots QC plot and spatial pattern of spots

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @param threshold lower limit on library size

#' @return spe_sub Spatial Experiment object with added col data after lib size thresholding

threshold_lib_size <- function(spe_sub, threshold) {

  # continuous variable containing total number of genes at each spot

  hist(colData(spe_sub)$sum, breaks = 20)

  # set filtering threshold for library size

  print(plotQC(spe_sub, type = "scatter", 

         metric_x = "cell_count", metric_y = "sum", 

         threshold_y = threshold))

  # check number of spots below threshold

  qc_lib_size <- colData(spe_sub)$sum < threshold

  print(table(qc_lib_size))

  # look for spatial pattern in discarded spots

  colData(spe_sub)$qc_lib_size <- qc_lib_size

  # plot spots

  threshold_row <- colData(spe_sub)[colData(spe_sub)$qc_lib_size,][,'array_row']

  threshold_col <- colData(spe_sub)[colData(spe_sub)$qc_lib_size,][,'array_col']

  print(ggplot(data.frame(threshold_row, threshold_col), aes(x=threshold_row, y=threshold_col)) + geom_point() + ggtitle("Discarded Spots (Low Library Size)"))

  spe_sub

}

#' Sets a lower filtering threshold for number of expressed genes, plots QC plot and spatial pattern of spots

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @param threshold lower limit on number of expressed genes

#' @return spe_sub Spatial Experiment object with added col data after expressed gene thresholding

threshold_expressed_genes <- function(spe_sub, threshold) {

  # continuous variable containing total number of genes at each spot

  hist(colData(spe_sub)$detected, breaks = 20)

  # set filtering threshold for library size at 500

  print(plotQC(spe_sub, type = "scatter", 

         metric_x = "cell_count", metric_y = "detected", 

         threshold_y = threshold))

  # check number of spots below threshold

  qc_detected <- colData(spe_sub)$detected < threshold

  print(table(qc_lib_size))

  # look for spatial pattern in discarded spots

  colData(spe_sub)$qc_detected <- qc_detected

  # plot spots

  threshold_row <- colData(spe_sub)[colData(spe_sub)$qc_detected,][,'array_row']

  threshold_col <- colData(spe_sub)[colData(spe_sub)$qc_detected,][,'array_col']

  print(ggplot(data.frame(threshold_row, threshold_col), aes(x=threshold_row, y=threshold_col)) + geom_point() + ggtitle("Discarded Spots (Low Expressed Genes)"))

  spe_sub

}

#' Sets an upper threshold on the proportion of mitochondrial reads

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @param threshold upper limit on proportional of mitochondrial reads

#' @return spe_sub Spatial Experiment object with added col data after mitochondrial read thresholding

threshold_mito_reads <- function(spe_sub, threshold) {

  # histogram of mitochondrial read proportions

  hist(colData(spe_sub)$subsets_mito_percent, breaks = 20)

  # plot mitochondrial read proportion vs. number of cells per spot

  print(plotQC(spe_sub, type = "scatter", 

         metric_x = "cell_count", metric_y = "subsets_mito_percent", 

         threshold_y = threshold))

  # select QC threshold for mitochondrial read proportion

  qc_mito <- colData(spe_sub)$subsets_mito_percent > threshold

  print(table(qc_mito))

  # look for spatial pattern in discarded spots

  colData(spe_sub)$qc_mito <- qc_mito

  # plot spots

  threshold_row <- colData(spe_sub)[colData(spe_sub)$qc_mito,][,'array_row']

  threshold_col <- colData(spe_sub)[colData(spe_sub)$qc_mito,][,'array_col']

  print(ggplot(data.frame(threshold_row, threshold_col), aes(x=threshold_row, y=threshold_col)) + geom_point() + ggtitle("Discarded Spots (High Mitochondrial Reads)"))

  spe_sub

}

#' Sets an upper threshold on the number of cells per spot

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @param threshold upper limit on number of cells per spot

#' @return spe_sub Spatial Experiment object with added col data after cell number thresholding

threshold_num_cells <- function(spe_sub, threshold) {

  # histogram of number of cells per spot

  hist(colData(spe_sub)$cell_count, breaks = 20)

  # distribution of cells per spot

  tbl_cells_per_spot <- table(colData(spe_sub)$cell_count)

  # plot number of expressed genes vs. number of cells per spot

  print(plotQC(spe_sub, type = "scatter", 

         metric_x = "cell_count", metric_y = "detected", 

         threshold_x = threshold))

  # select QC threshold for number of cells per spot

  qc_cell_count <- colData(spe_sub)$cell_count > threshold

  print(table(qc_cell_count))

  colData(spe_sub)$qc_cell_count <- qc_cell_count

  # plot spots

  threshold_row <- colData(spe_sub)[colData(spe_sub)$qc_cell_count,][,'array_row']

  threshold_col <- colData(spe_sub)[colData(spe_sub)$qc_cell_count,][,'array_col']

  print(ggplot(data.frame(threshold_row, threshold_col), aes(x=threshold_row, y=threshold_col)) + geom_point())

  spe_sub

}

#' Discard spots based on previous QC metrics

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @return spe_sub Spatial Experiment object with low quality spots removed

discard_spots <- function(spe_sub) {

  # get discarded spots from each QC operation

  qc_lib_size <- colData(spe_sub)$qc_lib_size

  qc_detected <- colData(spe_sub)$qc_detected

  qc_mito <- colData(spe_sub)$qc_mito

  qc_cell_count <- colData(spe_sub)$qc_cell_count

  # number of discarded spots for each metric

  apply(cbind(qc_lib_size, qc_detected, qc_mito, qc_cell_count), 2, sum)

  # combined set of discarded spots

  discard <- qc_lib_size | qc_detected | qc_mito | qc_cell_count

  print(table(discard))

  # store in object

  colData(spe_sub)$discard <- discard

  # check spatial pattern of combined set of discarded spots

  # plot spots

  threshold_row <- colData(spe_sub)[colData(spe_sub)$discard,][,'array_row']

  threshold_col <- colData(spe_sub)[colData(spe_sub)$discard,][,'array_col']

  print(ggplot(data.frame(threshold_row, threshold_col), aes(x=threshold_row, y=threshold_col)) + geom_point())

  spe_sub <- spe_sub[, !colData(spe_sub)$discard]

  print(dim(spe_sub))

  spe_sub

}

#' Discard uninformative genes from the Spatial Experiment object

#' @param spe_sub Spatial Experiment object filtered to only contain one sample

#' @param var_threshold lower threshold for amount of variance willing to tolerate among the spots for a gene

#' @return spe_sub Spatial Experiment object with uninformative genes removed

discard_uninformative_genes <- function(spe_sub, var_threshold=1e-4) {

  logcounts <- t(logcounts(spe_sub))

  mito <- (rowData(spe_sub)$gene_name %in% rowData(spe_sub)$gene_name[is_mito])

  zero_expression <- apply(logcounts, 2, function(x) length(unique(x)) == 1)

  low_var <- apply(logcounts, 2, function(x) var(x) < 1e-4) 

  discard <- zero_expression | low_var | mito

  spe_sub <- spe_sub[!discard,]

  spe_sub

}

```

### Driver Code

In the following block, data is first fetched from SpatialLIBD and filtered to only include a single sample. Then, the sample is plotted and QC metrics are added to the SpatialExperiment object.

Then, the following four QC operations are performed to get rid of low quality spots. Note that all spots are already over tissue.

1.  Set a lower threshold over library size. We want to remove spots with a low number of UMIs counts per spot.

2.  Set a lower threshold on the number of expressed features (i.e. genes) per spot.

3.  Set an upper threshold on the mumber of mitochondrial reads per spot. A high proportion of mitochondrial reads can indicate cell damage.

4.  Set an upper threshold on the number of cells per spot. Spots with high cell counts have a low number of expressed genes, meaning that experiments failed.

After performing QC operations, low quality spots are removed from the SpatialExperiment object.

```{r QC, message=FALSE}

# Fetch data from SpatialLIBD

spe <- fetch_data(type='spe')

# Filter spatial experiment object to only include sample ID

sample_id <- "151673"

spe_sub <- spe[, spe$sample_id == sample_id]

# Plot image and add QC metrics to spe object

img <- plot_image(spe_sub)

spe_sub <- add_QC_metrics(spe_sub)

# Perform QC operations

spe_sub <- threshold_lib_size(spe_sub, 800)

spe_sub <- threshold_expressed_genes(spe_sub, 550)

spe_sub <- threshold_mito_reads(spe_sub, 28)

spe_sub <- threshold_num_cells(spe_sub, 18)

# Discard low quality spots determined by QC operations

spe_sub <- discard_spots(spe_sub)

spe_sub <- discard_uninformative_genes(spe_sub)

```