cegs_2016.Rpres

```{r setup, include=FALSE}
opts_chunk$set(cache=TRUE, cache.lazy = FALSE)
```

Single-cell RNA-Seq Analysis
========================================================
author: Timothy Tickle and Brian Haas
css: cegs_2017.css
date: January 2017

Cut and Paste
===
class:small-code

Please follow along using the cut_and_paste.txt that comes with this repo.   


If code is cut off here in the presentation, it will be available in full in the cut_and_paste.txt document.

Logistics
===
class:small-code

- R allows methodology written by others to be imported.
  - Leverage other code.
  - Make your code available to others.  

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Load libraries
library(dplyr) # Dataframe manipulation
library(MAST) # For DE
library(Matrix) # Sparse matrices
library(useful) # Corner function
library(vioplot) # Violin pots
library(Seurat) # Single cell Analysis

source("./src/dotplot.R")
source("./src/makeNiceMAST.R")

load("./data/bipolar_raw.Robj")
```

Representing Sparse Matrices
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# 24904 x 44994 (1 minute)
bipolar.data <- Matrix(as.matrix(bipolar_dge), sparse=TRUE)
# Memory use as a normal matrix
object.size(bipolar_dge)

# Memory use as a sparse matrix
# ~23 X less
object.size(bipolar.data)
bipolar_dge <- NULL
```

Create a Seurat Object
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Expected raw counts (non-normalized data)
# Can give log transformed data but do not transform in setup method
bipolar.seurat.raw <- new("seurat", raw.data=bipolar.data)
```

What is in a Seurat Object?
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Display the internal pieces of the Seurat Object
slotNames(bipolar.seurat.raw)
```

What is in a Seurat Object?
===
class:small-code

- So far, raw sparse matrix.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
head(bipolar.seurat.raw@raw.data)
```

What is in a Seurat Object?
===
class:small-code

- So far, raw sparse matrix.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
head(bipolar.seurat.raw@data)
head(bipolar.seurat.raw@ident)
head(bipolar.seurat.raw@var.genes)
```

What is in a Seurat Object?
===
class:small-code

- __var.genes__ Variable genes across cells.
- __data.info__ Misc info including complexity (nGene).
- __cell.names__ Column (cell) names.
- __gene.names__ Row (gene) names.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
?seurat
```

What are Our Genes?
===
class:small-code

```{r,eval=TRUE, echo=TRUE, tidy=FALSE}
# Gene names (row names)
head(row.names(bipolar.seurat.raw@raw.data))
length(row.names(bipolar.seurat.raw@raw.data))
```

What are Our Cells?
===
class:small-code

```{r,eval=TRUE, echo=TRUE, tidy=FALSE}
# Column names # 24904 x 44994
# Sample / Cell names / Barcodes
head(colnames(bipolar.seurat.raw@raw.data),3)
length(colnames(bipolar.seurat.raw@raw.data))
dim(bipolar.seurat.raw@raw.data)
```

How to Show Counts?
===
class:small-code

```{r,eval=FALSE, echo=TRUE,tidy=FALSE}
# Only the corner
# The full data will be too large to see
corner(as.matrix(bipolar.seurat.raw@raw.data))
```

How to Show Counts?
===
class:small-code

```{r,eval=TRUE, echo=FALSE,tidy=FALSE}
# Only the corner
# The full data will be too large to see
corner(as.matrix(bipolar.seurat.raw@raw.data))
```

How Many Expressed Genes?
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Plot genes per cell
# How many genes expressed per cells
# 3 Minutes
complexity.per.cell <- apply(bipolar.seurat.raw@raw.data,
                             2, function(x) sum(x>0))
# Mean count per cell.
# 3 minutes
mean.count.per.cell <- apply(bipolar.seurat.raw@raw.data,
                             2, function(x) mean(x))
# Gene prevalence
# 5 minutes
gene.prevalence <- apply(bipolar.seurat.raw@raw.data,
                         1, function(x) sum(x>0))
```

Store Technical Batches for Later
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Get batches based on cell names
technical_batches <- sapply(colnames(bipolar.seurat.raw@raw.data),
                      FUN=function(x){strsplit(x,"_",fixed=TRUE)[[1]][1]})

# Turn to numbers and add cell names to them
technical_batches <- as.numeric(as.factor(technical_batches))
names(technical_batches) <- colnames(bipolar.seurat.raw@raw.data)
```

How Many Expressed Genes?
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Plot genes per cell
# How many genes expressed per cell
# Violin plot
vioplot(complexity.per.cell)

# Add points
stripchart(complexity.per.cell, add=TRUE, vertical=TRUE, method="jitter", jitter=0.3, pch=20, col="#00000033")

# Add lines for reference
abline(h=500, col="orange")
abline(h=3000, col="green")

# Add text
title("Study Complexity")
axis(side=1,at=1,labels=c("Study"))
```

How Many Expressed Genes?
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Plot genes per cell
# How many genes expressed per cell
vioplot(complexity.per.cell)
stripchart(complexity.per.cell, add=TRUE, vertical=TRUE, method="jitter", jitter=0.3, pch=20, col="#00000033")
abline(h=500, col="orange")
abline(h=3000, col="green")
title("Study Complexity")
axis(side=1,at=1,labels=c("Study"))
```

How Many Expressed Genes?
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Plot cell complexity, sorted by complexity
plot(1:length(complexity.per.cell),
     sort(complexity.per.cell),
     xlab="Cells ranked by complexity (Less to More)",
     ylab="Complexity")
# Add lines
abline(h=500, col="orange")
abline(h=3000, col="green")
title("Complexity from Less to More")
```

How Many Expressed Genes?
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
plot(1:length(complexity.per.cell), sort(complexity.per.cell))
abline(h=500, col="orange")
abline(h=3000, col="green")
```

How Many Expressed Genes By Batch?
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Plot genes per cell
# Complexity per technical batch
plot(complexity.per.cell ~ jitter(technical_batches,2))

# Add line
abline(h=500, col="orange")
abline(h=3000, col="green")

# Add text
title("Study Complexity")
axis(side=1,at=1,labels=c("Study"))
```

How Many Expressed Genes By Batch?
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Plot genes per cell
# How many genes expressed per cell
plot(complexity.per.cell ~ jitter(technical_batches,2))
abline(h=500, col="orange")
abline(h=3000, col="green")
title("Study Complexity")
axis(side=1,at=1,labels=c("Study"))
```

Identifying Outlier Cells?
===
class:small-code

- Cells that are unusually simple (or no counts).
- Cells that are unusually complex (doublets?).
- We will filter in a standard way, we are just describing the data.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Complexity by mean expression
plot(complexity.per.cell, mean.count.per.cell)

# Add lines
abline(v=500, col="orange")
abline(v=3000, col="green")
```

Identifying Outlier Cells?
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
plot(complexity.per.cell, mean.count.per.cell)
abline(v=500, col="orange")
abline(v=3000, col="green")
```

Identifying Noisey Genes?
===
class:small-code

- People tend to filter very conservatively.
- Remember this is in log space.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# How often genes are seen in cells
hist(log2(gene.prevalence), breaks=100)

# add line
abline(v=log(30), col="orange")
```

Identifying Noisey Genes?
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
hist(log2(gene.prevalence), breaks=100)
abline(v=log(30), col="orange")
```

Filter Cells: Removing the Outliers and Noise
===
class:small-code

- Genes must be in 30 cells.
- Cells must have alteast 500 genes.
- Scaled by 1000 (Total Sum Scaling).

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# From 24904 x 44994 to 14575 x 30673
# 5 minutes
bipolar.seurat <- Setup(bipolar.seurat.raw,
                    min.cells=30, min.genes=500,
                    do.logNormalize=TRUE,
                    total.expr=1e4,
                    project="Tutorial")
bipolar.seurat.raw <- NULL
```

Seurat: Updated Slots
===
class:small-code

```{r,echo=TRUE, eval=TRUE, tidy=FALSE}
dim(bipolar.seurat@raw.data)
```

Seurat: Updated Slots
===
class:small-code

```{r,echo=TRUE, eval=TRUE, tidy=FALSE}
dim(bipolar.seurat@data)
```

Seurat: Updated Slots
===
class:small-code

```{r,echo=TRUE, eval=TRUE, tidy=FALSE}
head(bipolar.seurat@data.info)
```

Seurat: Filtering on Metadata
===
class:small-code

- Filtering on mitochondrial reads in Seurat.
```{r, echo=TRUE, eval=TRUE, tidy=FALSE}
# Get gene names
mito.gene.names <- grep("^mt-", rownames(bipolar.seurat@data), value=TRUE)
```

Seurat: Filtering on Metadata
===
class:small-code

- Filtering on mitochondrial reads in Seurat.
```{r, echo=TRUE, eval=TRUE, tidy=FALSE}
# Get TSS normalized mitochodrial counts
# Cell total expression
col.total.counts <- Matrix::colSums(expm1(bipolar.seurat@data))

# Scale each count by the cell total
mito.percent.counts <- Matrix::colSums(expm1(bipolar.seurat@data[mito.gene.names, ]))/col.total.counts
```

Seurat: Filtering on Metadata
===
class:small-code

- Filtering on mitochondrial reads in Seurat.
```{r, echo=TRUE, eval=TRUE, tidy=FALSE}
# Add to seurat object as a metadata
bipolar.seurat <- AddMetaData(bipolar.seurat, mito.percent.counts, "percent.mitochodrial")

# Check
head(bipolar.seurat@data.info)
```

Seurat: Filtering on Metadata
===
class:small-code

- Plot current metadata in Seurat Object.
 - Number gene.
 - Number UMI.
 - Percent mitochondrial counts (batch affect).

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
VlnPlot(bipolar.seurat, "nGene")
VlnPlot(bipolar.seurat, "nUMI")
VlnPlot(bipolar.seurat, "percent.mitochodrial")
VlnPlot(bipolar.seurat, c("nGene", "nUMI", "percent.mitochodrial"), nCol=3)
```

Seurat: Filtering on Metadata
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "nGene")
```

Seurat: Filtering on Metadata
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "nUMI")
```

Seurat: Filtering on Metadata
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "percent.mitochodrial")
```

Seurat: Filtering on Metadata
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, c("nGene", "nUMI", "percent.mitochodrial"), nCol=3)
```

Seurat: Filtering on Metadata
===
class:small-code

- Outlier percent mitochondria are very low expression.
- Very high expression have low percent mitochondrial reads.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
GenePlot(bipolar.seurat, "nUMI", "percent.mitochodrial")
```

Seurat: Filtering on Metadata
===
class:small-code

- Outlier percent mitochondria are very low expression.
- Very high expression have low percent mitochondrial reads.

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
GenePlot(bipolar.seurat, "nUMI", "percent.mitochodrial")
```

Seurat: Filtering on Metadata
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Filter libraries with more than 10% mitochondrially derived transcripts
# From 14575 x 30673 to 14575 x 27489
dim(bipolar.seurat@data)
# 1 Minute
# Filter max complexity
bipolar.seurat <- SubsetData(bipolar.seurat, subset.name="nGene", accept.high=3000)

# Filter percent percent mitochondrial reads
bipolar.seurat <- SubsetData(bipolar.seurat, subset.name="percent.mitochodrial", accept.high=0.1)
```

Seurat: Filtering on Metadata
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
dim(bipolar.seurat@data)
table(bipolar.seurat@data.info[["orig.ident"]])
# Original standard matrix 8,971,430,520 bytes
object.size(bipolar.seurat)
```

Saving as an R Object
===
class:small-code

Saving the Seurat object.
- Contains all manipulation so far.
- Can be loaded or shared.
- Does not contain environment.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# 3.16 GB
# save(bipolar.seurat, file="./data/bipolar_seurat.Robj")

# load("./data/bipolar_seurat.Robj")
```

Saving as Text Files
===
class:small-code

You may need to export data to import into other applications.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Log-scale expression matrix
# 5 minutes (1.2 GB)
# write.table(as.matrix(bipolar.seurat@data), file="bipolar_data.txt")

# Study metadata (1.2 MB)
# write.table(bipolar.seurat@data.info, file="bipolar_metadata.txt")
```

Seurat: Viewing Specific Genes
===
class:small-code

Plotting genes across groups.
- Check the identity of the cells!!!

```{r,eval=FALSE, echo=TRUE, tidy=FALSE}
VlnPlot(bipolar.seurat, "Prkca")
VlnPlot(bipolar.seurat, "Prkca", group.by="orig.ident")
```

Seurat: Viewing Specific Genes
===
class:small-code

- Notice many zeros.

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "Prkca")
```

Seurat: Viewing Specific Genes
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "Prkca", group.by="orig.ident")
```

Seurat: Plotting Genes vs Genes
===
class:small-code

```{r,eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
# Plot a gene vs a gene
GenePlot(bipolar.seurat,"Isl1","Grm6",cex.use=1)
```

Seurat: Plotting Genes vs Genes
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Plot a gene vs a gene
GenePlot(bipolar.seurat,"Isl1","Grm6",cex.use=1)
```

Seurat: Batch Affect Correction
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Create batch affect (5 and 6 as a group)
batch.effect <- technical_batches[row.names(bipolar.seurat@data.info)]
batch.effect[batch.effect %in% c(1,2,3,4)] <- 1
batch.effect[batch.effect %in% c(5,6)] <- 2

# Add batch affect
bipolar.seurat <- AddMetaData(bipolar.seurat, batch.effect, "batch.effect")

# Regress on metadata
bipolar.seurat <- RegressOut(bipolar.seurat, latent.vars="batch.effect")
```

Seurat: Select Variable Genes
===
class:small-code

- We are going to focus on highly variable genes.
  - Average expression (X) and dispersion (SD)(Y).
  - Bins genes (20).
  - Z-score for dispersion.
- fxn.x and fxn.y allow one to change the measurements.
  - Cut.offs are high and low, X and Y.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# 2 minutes
bipolar.seurat <- MeanVarPlot(bipolar.seurat,
                              fxn.x=expMean,
                              fxn.y=logVarDivMean,
                              x.low.cutoff=0.0125,
                              x.high.cutoff=3,
                              y.cutoff=0.5, do.contour=FALSE,
                              do.plot=FALSE)
length(bipolar.seurat@var.genes)
```

Seurat: Performing PCA
===
class:small-code

- Calculating PCA with the highly variable genes.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE, fig.align="center"}
# 9 minutes
bipolar.seurat <- PCA(bipolar.seurat,
                      pc.genes=bipolar.seurat@var.genes,
                      do.print=FALSE)
# save(bipolar.seurat, file="./data/bipolar_seurat_pca.Robj")
```

Seurat: Performing PCA
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE, fig.align="center"}
# Calculate PCA projection
bipolar.seurat <- ProjectPCA(bipolar.seurat)
```

Seurat: Performing PCA
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE, fig.align="center"}
# Can plot top genes for top components
PrintPCA(bipolar.seurat,
         pcs.print=1:2,
         genes.print=5,
         use.full=TRUE)
```

Seurat: PCA Visualizations
===
class:small-code

- Top 30 genes associated with the first two components.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
VizPCA(bipolar.seurat, pcs.use=1:2)
```

Seurat: PCA Visualizations
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VizPCA(bipolar.seurat, pcs.use=1:2)
```

Seurat: Check your Batch Correction
===
class:small-code

How can we use component loadings?
- Can check components to see if genes show up that are enriched in components.
- Using the scores you can perform enrichment analysis to decribe the component.

Seurat: PCA Visualizations
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
PCAPlot(bipolar.seurat, 1, 2)
```

Seurat: PCA Visualizations
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE}
PCAPlot(bipolar.seurat, 1, 2)
```

Seurat: PCA Visualizations
===
class:small-code

- Plot top 30 genes in top 100 cells for PC1.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
PCHeatmap(bipolar.seurat, pc.use=1, cells.use=100, do.balanced=TRUE)
```

Seurat: PCA Visualizations
===
class:small-code

- Plot top 30 genes in top 100 cells for PC1.

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
PCHeatmap(bipolar.seurat, pc.use=1, cells.use=100, do.balanced=TRUE)
```

Seurat: Choosing Components
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Time Intensive
# Jackstraw
# bipolar.seurat <- JackStraw(bipolar.seurat, num.replicate = 100, do.print = FALSE)
```

Seurat: Choosing Components
===
class:small-code

How do we choose how many components to use?
 - When there is diminishing returns to include it.
 - Selection is performed more liberally in our setting.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Time Efficient but more ad hoc
# Scree (elbow) plot
# 37 selected
PCElbowPlot(bipolar.seurat, num.pc=40)
```

Seurat: Choosing Components
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
PCElbowPlot(bipolar.seurat, num.pc=40)
```

Seurat: Run t-SNE
===
class:small-code

- Calculate and plot t-SNE on PC 1- 37.
- Can use Barnes-hut implementation.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Calculate t-SNE Ordination
# 7 minutes
bipolar.seurat <- RunTSNE(bipolar.seurat,
                          dims.use=1:37,
                          do.fast=TRUE)
```

Seurat: Plot t-SNE
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE}
# Plot
TSNEPlot(bipolar.seurat)
```

Seurat: Plot t-SNE
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Plot
TSNEPlot(bipolar.seurat)
```

Seurat: Plot t-SNE
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Plot
TSNEPlot(bipolar.seurat)
```

---

This is not PCA
- Distance is best understood in close neighbors.
- The measurement of distance is difficult to understand.

Seurat: Side by side
===
class:small-code

**PCA**
```{r, eval=TRUE, echo=FALSE, tidy=FALSE}
PCAPlot(bipolar.seurat, 1, 2)
```

---

**t-SNE**
```{r, eval=TRUE, echo=FALSE, tidy=FALSE}
TSNEPlot(bipolar.seurat)
```

Seurat: Calculate Clusters
===
class:small-code

- Determine subclusters for the plot.
- Separate graph based approach.
   - is not aware of the t-SNE projection.
- Using PC 1-37
```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# 14 minutes
# bipolar.seurat <- FindClusters(bipolar.seurat,
#                            pc.use=1:37,
#                            resolution=0.6,
#                            print.output=0,
#                            save.SNN=TRUE)
# save(bipolar.seurat, file="./data/bipolar_cluster_seurat.Robj")
load("./data/bipolar_cluster_seurat.Robj")
```

Seurat: Viewing de novo Groupings
===
class:small-code

```{r,eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
TSNEPlot(bipolar.seurat)
```

Seurat: Viewing de novo Groupings
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
TSNEPlot(bipolar.seurat)
```

Seurat: Viewing de novo Groupings
===
class:small-code

```{r,eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
TSNEPlot(bipolar.seurat, do.label=TRUE)
```

Seurat: Viewing de novo Groupings
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
TSNEPlot(bipolar.seurat, do.label=TRUE)
```

Seurat: Plotting Genes Through Clusters
===
class:small-code

Now that we have subclusters of cell populations, plotting genes through subclusters is identical to before.
- Seurat stores and groups by subclusters automatically.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "Prkca")
```

Seurat: Plotting Genes Through Clusters
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "Prkca")
```

Seurat: Plotting Genes Through Clusters
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "Glul")
```

Seurat: Plotting Genes Through Clusters
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
VlnPlot(bipolar.seurat, "Glul")
```

Seurat: Plotting Genes on Clusters
===
class:small-code

You can also gene expression through out the cell ordination.
- Marker genes can help identify cell groups.
  - Rod Bipolar Cells (Prkca+ Scgn-).
  - Muller Glia (Glul+ Prkca- Scgn-).
  - Cone Bipolar Cells (Prkca- Scgn+).
  - On Cone Bipolar Cells (Prkca- Scgn+ Grm6+).
- Metadata can help visualize batch affects.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
FeaturePlot(bipolar.seurat,
            c("Prkca","Glul", "Scgn", "Grm6"),
            cols.use = c("grey","blue"))
```

Seurat: Plotting Genes on Clusters
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
FeaturePlot(bipolar.seurat,
            c("Prkca","Glul", "Scgn", "Grm6"),
            cols.use = c("grey","blue"))
```

QC the Clusters!
===
class:small-code

```{r,eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
FeaturePlot(bipolar.seurat,
            "nGene",
            cols.use = c("grey","blue"))
FeaturePlot(bipolar.seurat,
            "percent.mitochodrial",
            cols.use = c("grey","blue"))
```

QC the Clusters!
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
FeaturePlot(bipolar.seurat,
            "nGene",
            cols.use = c("grey","blue"))
```

QC the Clusters!
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
FeaturePlot(bipolar.seurat,
            "percent.mitochodrial",
            cols.use = c("grey","blue"))
```

QC the Clusters!
===
class:small-code

Check for your batch affect.
- We are going to make a fake batch affect (site) and plot this as an example of how one can visualize unwanted signal.

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
# Making Fake Data
fake.sites <- as.integer(bipolar.seurat@ident %in% c(0))
names(fake.sites) <- names(bipolar.seurat@ident)

# Add metadata
bipolar.seurat <- AddMetaData(bipolar.seurat, fake.sites, "site")

# Plot feature
FeaturePlot(bipolar.seurat, c("site"), cols.use = c("green","orange"))
```

QC the Clusters!
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Making Fake Data
fake.sites <- as.integer(bipolar.seurat@ident %in% c(0))
names(fake.sites) <- names(bipolar.seurat@ident)

# Add metadata
bipolar.seurat <- AddMetaData(bipolar.seurat, fake.sites, "site")

# Plot feature
FeaturePlot(bipolar.seurat, c("site"), cols.use = c("green","orange"))
```

Seurat: Getting your labels
===
class:small-code

```{r,eval=FALSE, echo=TRUE, tidy=FALSE}
cell.labels <- bipolar.seurat@ident
corner(cell.labels)
```

Seurat: Getting your labels
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE}
cell.labels <- bipolar.seurat@ident
corner(cell.labels)
```

Seurat: Differential Expression
===
class:small-code

- Default is one cluster against many tests.
  - Can specify an ident.2 test between clusters.
- Adding speed by excluding tests.
  - Min.pct - controls for sparsity.
    - Min percentage in a group.
  - Thresh.test - must have this difference in averages.

```{r,eval=FALSE, echo=TRUE, tidy=FALSE}
# 2 minutes
cluster5.markers <- FindMarkers(bipolar.seurat, ident.1=5, min.pct = 0.25)
head(cluster5.markers, 30)
```

Seurat: Differential Expression
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE}
# 2 minutes
cluster5.markers <- FindMarkers(bipolar.seurat, ident.1=5, min.pct = 0.25)
head(cluster5.markers, 30)
```

Seurat: Differential Expression
===
class:small-code

- Find all cluster against all others.

```{r,eval=FALSE, echo=TRUE, tidy=FALSE}
# bipolar.markers <- FindAllMarkers(bipolar.seurat, only.pos=TRUE, min.pct=0.25, thresh.use=0.25)
# bipolar.markers %>% group_by(cluster) %>% top_n(2, avg_diff)
```

Seurat: DE Tests
===
class:small-code

**bimod:** Tests differences in mean and proportions.   

**roc:** Uses AUC like definition of separation.   

**t:** Student's T-test.   

**tobit:** Tobit regression on a smoothed data.   

Seurat: Plotting DE Genes
===
class:small-code

```{r,eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
plot.genes <- head(rownames(cluster5.markers),30)
DoHeatmap(bipolar.seurat,
          genes.use=plot.genes,
          order.by.ident=TRUE,
          slim.col.label=TRUE,
          remove.key=TRUE)
```

Seurat: Plotting DE Genes
===
class:small-code

```{r,eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
plot.genes <- head(rownames(cluster5.markers),30)
DoHeatmap(bipolar.seurat,
          genes.use=plot.genes,
          order.by.ident=TRUE,
          slim.col.label=TRUE,
          remove.key=TRUE)
```

Dot Plots: Plotting Genes Through Clusters
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
dot.plot(bipolar.seurat,
         features.use=plot.genes)
```

Dot Plots: Plotting Genes Through Clusters
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
dot.plot(bipolar.seurat,
         features.use=plot.genes)
```

Dot Plots: Replicating Bipolar Results
===
class:small-code

```{r, eval=FALSE, echo=TRUE, tidy=FALSE, fig.align="center"}
# Genes in tutorial
plot.genes = c("Vsx2","Otx2","Scgn","Isl1","Grm6",
               "Apoe","Pax6","Rho","Arr3","Tacr3",
               "Syt2","Neto1","Irx6","Prkar2b",
               "Grik1","Kcng4","Cabp5","Vsx1","Prkca")

# Cell groups shown
only.use.groups = c(7,9,10,12,8,14,3,13,6,11,5,4,16,0)

# Plot
dot.plot(bipolar.seurat,
         features.use=plot.genes,
         subset.group=only.use.groups)
```

Dot Plots: Replicating Bipolar Results
===
class:small-code

```{r, eval=TRUE, echo=FALSE, tidy=FALSE, fig.align="center"}
# Genes in tutorial
plot.genes = c("Vsx2","Otx2","Scgn","Isl1","Grm6",
               "Apoe","Pax6","Rho","Arr3","Tacr3",
               "Syt2","Neto1","Irx6","Prkar2b",
               "Grik1","Kcng4","Cabp5","Vsx1","Prkca")

# Cell groups shown
only.use.groups = c(7,9,10,12,8,14,3,13,6,11,5,4,16,0)

# Plot
dot.plot(bipolar.seurat,
         features.use=plot.genes,
         subset.group=only.use.groups)
```

MAST: Creating a Test Comparison
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Make a covariate to inform the test
test.5 <- as.character(bipolar.seurat@ident)
test.5[test.5 != "5"] <- "wild"
test.5[test.5 == "5"] <- "test"
names(test.5) <- names(bipolar.seurat@ident)
test.5 <- factor(test.5, levels=c("wild","test"))
```

MAST: Prepping the Data
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Load into Mast Object
# Mast object
# Data frame of cell-level covariates
# Data frame of feature-level covariates
bipolar.mast <- FromMatrix(as.matrix(bipolar.seurat@data),
                  cbind(bipolar.seurat@data.info,test.5),
                  data.frame(primerid=rownames(bipolar.seurat@data)))
```

MAST: Create the Model
===
class:small-code

- Create the model, here we load data for time.
- Perform test for our contrast.
- Format results for easier use.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
# Run test
# 1 hr
# zlm.test.5 <- zlm.SingleCellAssay(~ test.5, bipolar.mast)
load("./data/mast_zlm.Robj")

df.test.5 <- data.frame(summary(zlm.test.5, doLRT="test.5test"))
summary.mast <- makeNice(df.test.5, val="test.5test")
```

MAST: View Results
===
class:small-code

- Results from test.
- Some overlapping genes.

```{r, eval=TRUE, echo=TRUE, tidy=FALSE}
summary.mast[abs(summary.mast[,5])>=1 & !is.na(summary.mast[,5]),]
```

MAST: Plot Results
===
class:small-code

```{r, eval=TRUE, echo=TRUE, tidy=FALSE, fig.align="center"}
plot.genes = c("Scg2","Scgn","Cck","Vsx1","Guk1",
               "Gngt2","Pcp4l1","Pcp4","Pcp2")
dot.plot(bipolar.seurat,
         features.use=plot.genes)
```