QC.Rmd

---
title: "QC"
author: "Harvard Chan Bioinformatics Core"
date: "`r Sys.Date()`"
output:
   html_document:
      code_folding: hide
      df_print: paged
      highlights: pygments
      number_sections: false
      self_contained: true
      theme: default
      toc: true
      toc_float:
         collapsed: true
         smooth_scroll: true
params:
  project_file: ./reports/information.R
  seurat_fn: ./data/CTCL_PAT8_DEV3_LVL7.RDS
  results_dir: ./results
  min_transcripts: 100
  min_genes: 20
  min_novelty: .7
  umap_dim: approximateumap_8c6f278e.b9f4.4535.aeca.8955c1dff614_1
---

```{r, cache = FALSE, message = FALSE, warning=FALSE}
# This set up the working directory to this file so all files can be found
library(rstudioapi)
setwd(fs::path_dir(getSourceEditorContext()$path))
# NOTE: This code will check version, this is our recommendation, it may work
#.      other versions
stopifnot(R.version$major>= 4) # requires R4
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1") 
stopifnot(compareVersion(as.character(BiocManager::version()), "3.18")>=0)
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.0.0")>=0)
```

This code is in this ![](https://img.shields.io/badge/status-draft-grey) revision.

```{r load_params, echo = F}
source(params$project_file)
```

```{r load_libraries, cache = FALSE, message = FALSE, warning=FALSE, echo=FALSE,}

library(tidyverse)
library(Seurat)
library(bcbioR)
library(ggprism)
library(knitr)
library(tools)
library(qs)

colors=cb_friendly_cols(1:15)
ggplot2::theme_set(theme_prism(base_size = 14))
opts_chunk[["set"]](
    cache = F,
    cache.lazy = FALSE,
    dev = c("png", "pdf"),
    error = TRUE,
    highlight = TRUE,
    message = FALSE,
    prompt = FALSE,
    tidy = FALSE,
    warning = FALSE,
    echo = T, 
    fig.height = 4)

# set seed for reproducibility
set.seed(1234567890L)
```

```{r sanitize_datatable}
sanitize_datatable = function(df, ...) {
 # remove dashes which cause wrapping
 DT::datatable(df, ..., rownames=gsub("-", "_", rownames(df)),
                   colnames=gsub("-", "_", colnames(df)),
               filter = 'top')
}
```

# Overview

-   Project: `r project`
-   PI: `r PI`
-   Analyst: `r analyst`
-   Experiment: `r experiment`
-   Aim: `r aim`

```{r read rds}

# resave RDS as QS for faster reading/writing
seurat_qs_fn <- paste0(file_path_sans_ext(params$seurat_fn), '.qs')

if (!file.exists(seurat_qs_fn)){
  seurat <- readRDS(params$seurat_fn)
  qsave(seurat, seurat_qs_fn, preset = 'fast')
} else {
  seurat <- qread(seurat_qs_fn)
}

centroids <- data.frame(x = seurat$x_slide_mm, y = seurat$y_slide_mm, cell = seurat$cell_id)
cents <- CreateCentroids(centroids)
coords <- CreateFOV(coords = list(centroids = cents), type = "centroids")

```

```{r plot tissue}

ggplot(seurat@meta.data, aes(y = x_slide_mm, x = y_slide_mm)) + 
  geom_point(alpha = 0.05, size = 0.01) + facet_wrap(~Run_Tissue_name) + coord_equal() + 
  labs(title = "Cell coordinates in XY space")


```

# QC Plots {.tabset}

## nGenes

```{r plot n_genes hist}
ggplot(data = seurat[[]], aes(x = nFeature_RNA)) + geom_histogram(binwidth = 50) +
     geom_vline(xintercept = params$min_genes, col = "red", linetype = "dashed") +
     xlab("nGenes")
```

## nUMIs

```{r plot n_umis hist}
ggplot(data = seurat[[]], aes(x = nCount_RNA)) + geom_histogram(binwidth = 50) +
     geom_vline(xintercept = params$min_transcripts, col = "red", linetype = "dashed") +
     xlab("nUMIs") 
```

## nUMIs ranked
```{r plot n_umis ranked}
nUMIs_df <- seurat[[]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var = "barcode") %>%
  dplyr::select(c(barcode, nCount_RNA)) %>%
  dplyr::arrange(-nCount_RNA)

nUMIs_df$ranked_bc <- as.integer(rownames(nUMIs_df))

```

## novelty
```{r plot novelty}
seurat$novelty <- log10(seurat@meta.data$nFeature_RNA) / log10(seurat@meta.data$nCount_RNA)

novelty_df <- seurat[[]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var = "barcode") %>%
  dplyr::select(c(barcode, novelty)) %>%
  dplyr::arrange(-novelty)

novelty_df$ranked_bc <- as.integer(rownames(novelty_df))

ggplot() + geom_line(data = novelty_df, aes(x = ranked_bc, y = novelty), col = "red") +
  scale_y_continuous(trans = "log10") +
  geom_hline(yintercept = params$min_novelty,  linetype = "dashed") +
  ylab("novelty") + xlab("barcode")

```

# Filtering Low-Quality Cells

We discard cells that have less than `r params$min_genes` features and genes present in less than 3 cells. Additionally, we apply the following AtoMx QC flags to select/filter cells:

qcFlagsCellComplex = Pass. Filtering for complexity (nCount_RNA / nFeature_RNA) qcFlagsCellArea = Pass. Cell areas flagged as outliers by Grubbs test


```{r qc filtering}
counts <- LayerData(seurat, layer = "counts", assay = "RNA")

seurat_filtered <- CreateSeuratObject(counts = counts, meta.data = seurat@meta.data, min.cells = 3,
     min.features = params$min_genes)

selected_cells <- seurat_filtered[[]] %>%
     as.data.frame() %>%
     dplyr::filter(qcFlagsCellComplex == "Pass", 
                   qcFlagsCellArea == "Pass",
                   novelty > params$min_novelty
                   # qcFlagsFOV == "Pass"
                   ) %>%
     pull(cell_id)

seurat_filtered <- subset(seurat_filtered, cells = selected_cells)

n_before <- nrow(seurat[[]])
n_after <- nrow(seurat_filtered[[]])
```

There were `r n_before` cells before filtering and `r n_after` afterwards, for a total of
`r n_after/n_before * 100`% remaining

```{r plot pre/post qc}

pre <- ImageDimPlot(seurat) + NoLegend() + ggtitle("Pre-Filtering")
seurat_filtered[["FOV"]] <- subset(coords, cell = Cells(seurat_filtered))
post <- ImageDimPlot(seurat_filtered) + NoLegend() + ggtitle("Post-Filtering")

discarded_cells <- colnames(seurat)[!colnames(seurat) %in% colnames(seurat_filtered)]
seurat$selected <- seurat[[]] %>%
  as.data.frame() %>%
  dplyr::mutate(
    selected = case_when(
      cell_id %in% discarded_cells ~ "discarded", 
      TRUE ~ "selected")
    ) %>%
  pull(selected)

seurat$selected <- factor(seurat$selected, levels = c("selected", "discarded"))
discarded <- ImageDimPlot(seurat, group.by = "selected") + NoLegend() + ggtitle("Blue - discarded")
print(pre + post + discarded)

```

# Processing

```{r processing}

# perform processing (one time). if already done previously, load from file
processed_seurat_fn <- paste0(file_path_sans_ext(params$seurat_fn), '_processed.qs')
if (!file.exists(processed_seurat_fn)) {
  seurat_filtered <- SCTransform(seurat_filtered, assay = "RNA", clip.range = c(-10,10), verbose = FALSE)
  seurat_filtered <- NormalizeData(seurat_filtered, assay = "RNA")
  seurat_filtered <- RunPCA(seurat_filtered)
  seurat_filtered <- FindNeighbors(seurat_filtered, dims = 1:30)
  seurat_filtered <- RunUMAP(seurat_filtered, dims = 1:30)
  seurat_filtered <- FindClusters(seurat_filtered, resolution = 0.1, verbose = FALSE)
  qsave(seurat_filtered, processed_seurat_fn, preset = 'fast')
} else {
  seurat_filtered <- qread(processed_seurat_fn)
}

```

```{r plot umap before}

# TODO find colname of pre-filtering umap data in seurat object, use as params$umap_dim at top of file
DimPlot(seurat, reduction = paste(params$umap_dim),
     pt.size = 0.6) + labs(x = "umap_1", y = "umap_2", title = "Pre-Filtering")

```

```{r plot umap after}
DimPlot(seurat_filtered, reduction = "umap", pt.size = 0.6) +
     ggtitle("Post-Filtering")

```

```{r plot image clusters}
ImageDimPlot(seurat_filtered, axes = TRUE, crop = TRUE, combine = TRUE)

```

# Markers

## Cell Type Markers {.tabset}

```{r markers of interest, results = 'asis'}

## TODO replace with markers relevant to your project
markers_of_interest <- c('CD4', 'CD8A', 'CD8B', 'CD63', 'CD69', 'HBB')

for (marker in markers_of_interest) {
  cat("### ", marker, "\n")
  
  FeaturePlot(seurat_filtered, features = marker, max.cutoff = "q95", min.cutoff = "q05",
              reduction = "umap", pt.size = 0.6, order = T) %>%
    print()
 
  p <- ImageFeaturePlot(seurat_filtered, features = marker, max.cutoff = "q95",
                        size = 1, crop = TRUE, combine = FALSE)
  print(p)
  cat('\n')
}
```

# Cell Type Identification

## Azimuth

```{r}
# perform Azimuth cell type identification (once). if already done previously, load from file
azimuth_seurat_fn <- paste0(file_path_sans_ext(params$seurat_fn), '_azimuth.qs')
if (!file.exists(azimuth_seurat_fn)) {
  seurat_filtered_ann_pbmc <- RunAzimuth(seurat_filtered, assay = "RNA", reference = "pbmcref")
  qsave(seurat_filtered_ann_pbmc, azimuth_seurat_fn, preset = 'fast')

} else {
  seurat_filtered_ann_pbmc <- qread(azimuth_seurat_fn)
}
DimPlot(seurat_filtered_ann_pbmc, group.by = 'predicted.celltype.l1')
ImageDimPlot(seurat_filtered_ann_pbmc, axes = TRUE, crop = TRUE, combine = TRUE, group.by = 'predicted.celltype.l1')

```

## Seurat and GPT4
```{r}
markers <- FindAllMarkers(seurat_filtered)

# TODO: uncomment this chunk to get markers for pasting into GPT4
# markers %>%
#   dplyr::filter(avg_log2FC > 0) %>%
#   select(c("cluster", "gene")) %>%
#   group_by(cluster) %>%
#   slice(1:20) %>%
#   summarise(gene = paste(gene, collapse = ", "))

markers %>% sanitize_datatable()

## TODO replace with cell types identified by GPT4 for your markers
cluster_ids <- data.frame(cluster = as.factor(c(0:6)),
                          cell_type = as.factor(c('B-cells', 'Fibroblasts', 'Monocytes',
                                        'Cytotoxic T-cells', 'Endothelial Cells',
                                        'T-cells', 'Keratinocytes')))


seurat_filtered@meta.data <- left_join(seurat_filtered@meta.data, cluster_ids, by = c('seurat_clusters' = 'cluster'))
rownames(seurat_filtered@meta.data) <- seurat_filtered@meta.data$cell_id

ImageDimPlot(seurat_filtered, group.by = 'cell_type', axes = TRUE, crop = TRUE)


```