---
title: "Using sctransform in Seurat"
author: "Christoph Hafemeister & Rahul Satija"
date: '`r Sys.Date()`'
output: 
  html_document:
    highlight: pygments
---

```{r setup, include = FALSE}
library('Matrix')
library('ggplot2')
library('reshape2')
library('sctransform')
library('knitr')
knit_hooks$set(optipng = hook_optipng)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  digits = 2,
  tidy = TRUE,
  tidy.opts = list(width.cutoff=80),
  optipng = '-o 5 -strip all -quiet',
  fig.width=13, fig.height=6.5, dpi=100, out.width = '95%'
)
library('Seurat')
library('patchwork')
#old_theme <- theme_set(theme_classic(base_size=8))
# some of the vst steps can use multiple cores
# We use the Future API for parallel processing; set parameters here
future::plan(strategy = 'multicore', workers = 3)
options(future.globals.maxSize = 8 * 1024 ^ 3)
options(future.fork.enable = TRUE)
```

This vignette shows how to use the sctransform wrapper in Seurat.

Make sure the Seurat package is loaded

```{r eval=FALSE}
library(Seurat)
```

Load data and create Seurat object

```{r load_data, warning=FALSE, message=FALSE, cache = T}
pbmc_data <- Read10X(data.dir = "~/Downloads/pbmc3k_filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc_data)
```

For reference, we first apply the standard Seurat workflow, with log-normalization

```{r logtrans, warning=FALSE, message=FALSE, cache = T}
pbmc_logtransform <- pbmc
pbmc_logtransform <- NormalizeData(pbmc_logtransform, verbose = FALSE)
pbmc_logtransform <- FindVariableFeatures(pbmc_logtransform, verbose = FALSE) 
pbmc_logtransform <- ScaleData(pbmc_logtransform, verbose = FALSE) 
pbmc_logtransform <- RunPCA(pbmc_logtransform, verbose = FALSE) 
pbmc_logtransform <- RunUMAP(pbmc_logtransform,dims = 1:20, verbose = FALSE)
```

For comparison, we now apply sctransform normalization

```{r apply_sct, warning=FALSE, message=FALSE, cache = T}
# Note that this single command replaces NormalizeData, ScaleData, and FindVariableFeatures.
# Transformed data will be available in the SCT assay, which is set as the default after running sctransform
pbmc <- SCTransform(object = pbmc, verbose = FALSE)
```

Perform dimensionality reduction by PCA and UMAP embedding

```{r pca, cache = T}
# These are now standard steps in the Seurat workflow for visualization and clustering
pbmc <- RunPCA(object = pbmc, verbose = FALSE)
pbmc <- RunUMAP(object = pbmc, dims = 1:20, verbose = FALSE)
pbmc <- FindNeighbors(object = pbmc, dims = 1:20, verbose = FALSE)
pbmc <- FindClusters(object = pbmc, verbose = FALSE)
```

Visualize the clustering results on the sctransform and log-normalized embeddings.
```{r viz, cache = T}
pbmc_logtransform$clusterID <- Idents(pbmc)
Idents(pbmc_logtransform) <- 'clusterID'
plot1 <- DimPlot(object = pbmc, label = TRUE) + NoLegend() + ggtitle('sctransform') 
plot2 <- DimPlot(object = pbmc_logtransform, label = TRUE)
plot2 <- plot2 + NoLegend() + ggtitle('Log-normalization') 
plot1 + plot2
```

Users can individually annotate clusters based on canonical markers. However, the sctransform normalization reveals sharper biological distinctions compared to the log-normalized analysis. For example, note how clusters 0, 1, 9, and 11 (all distinct clusters of T cells), are blended together in log-normalized analyses. The sctransform analysis reveals:


* Clear separation of three CD8 T cell clusters (naive, memory, effector), based on CD8A, GZMK, CCL5 expression 
* Clear separation of three CD4 T cell clusters (naive, memory, IFN-activated) based on S100A4, CCR7, IL32, and ISG15 
* Additional developmental sub-structure in B cell cluster, based on TCL1A, FCER2
* Additional sub-structure within NK cell cluster (CD56dim vs. bright), based on XCL1 and FCGR3A

Visualize canonical marker genes as violin plots.

```{r fplot, cache = T}
VlnPlot(object = pbmc, features = c("CD8A","GZMK","CCL5","S100A4","S100A4","CCR7","ISG15","CD3D"), pt.size = 0.2,ncol = 4)
```

Visualize canonical marker genes on the sctransform embedding.

```{r fplot2, cache = T}
FeaturePlot(object = pbmc, features = c("CD8A","GZMK","CCL5","S100A4","S100A4","CCR7"), pt.size = 0.2,ncol = 3)
```

```{r fplot3, cache = T}
FeaturePlot(object = pbmc, features = c("CD3D","ISG15","TCL1A","FCER2","XCL1","FCGR3A"), pt.size = 0.2,ncol = 3)
```
