---
title: "Detecting hidden heterogeneity in single cell RNA-Seq data"
author: "Donghyung Lee"
date: "May 10th, 2018"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{"Detecting hidden heterogeneity in single cell RNA-Seq data"}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


The iasva package can be used to detect hidden heterogenity within bulk or 
single cell sequencing data. To illustrate how to use the iasva package for
heterogenity detection, we use real-world single cell RNA sequencing (scRNA-Seq)
data obtained from human pancreatic islet samples 
([Lawlor et. al., 2016](http://genome.cshlp.org/content/early/2017/01/16/gr.212720.116)). 

## Load packages
```{r load_packages, echo=TRUE, message=FALSE}
library(irlba) 
library(iasva)
library(sva)
library(Rtsne)
library(pheatmap)
library(corrplot)
library(DescTools)
library(RColorBrewer)
library(SummarizedExperiment)
set.seed(100)
color.vec <- brewer.pal(3, "Set1")
```

## Load the islet single cell RNA-Seq data
To illustrate how IA-SVA can be used to detect hidden heterogeneity within
a homogenous cell population (i.e., alpha cells), we use read counts of alpha
cells from healthy (non-diabetic) subjects (n = 101).   
```{r load_data, echo=TRUE}
counts_file <- system.file("extdata", "iasva_counts_test.Rds",
                         package = "iasva")
# matrix of read counts where genes are rows, samples are columns
counts <- readRDS(counts_file)
# matrix of sample annotations/metadata
anns_file <- system.file("extdata", "iasva_anns_test.Rds",
                         package = "iasva")
anns <- readRDS(anns_file)

```

## Calculate geometric library size, i.e., library size of log-transfromed read counts. 
It is well known that the geometric library size
(i.e., library size of log-transfromed read counts) or proportion of expressed
genes in each cell explains a very large portion of variability of scRNA-Seq 
data ([Hicks et. al. 2015 BioRxiv](http://biorxiv.org/content/early/2015/08/25/025528),
[McDavid et. al. 2016 Nature Biotechnology](http://www.nature.com/nbt/journal/v34/n6/full/nbt.3498.html)). 
Frequently, the first principal component of log-transformed scRNA-Seq read
counts is highly correlated with the geometric library size (r ~ 0.9). 
Here, we calculate the geometric library size vector, which will be used
as a known factor in the IA-SVA algorithm.
```{r geno_lib_size, echo=TRUE, fig.width=7, fig.height=4}
geo_lib_size <- colSums(log(counts + 1))
barplot(geo_lib_size, xlab = "Cell", ylab = "Geometric Lib Size", las = 2)
lcounts <- log(counts + 1)

# PC1 and Geometric library size correlation
pc1 <- irlba(lcounts - rowMeans(lcounts), 1)$v[, 1]
cor(geo_lib_size, pc1)
```

## Run IA-SVA
Here, we run IA-SVA using patient_id and geo_lib_size as known factors
and identify five hidden factors. SVs are plotted in a pairwise fashion to
uncover which SVs can seperate cell types. 
```{r run_iasva, echo=TRUE, fig.width= 7, fig.height=6}
set.seed(100)
patient_id <- anns$Patient_ID
mod <- model.matrix(~patient_id + geo_lib_size)
# create a summarizedexperiment class
summ_exp <- SummarizedExperiment(assays = counts)
iasva.res<- iasva(summ_exp, mod[, -1],verbose = FALSE, 
                  permute = FALSE, num.sv = 5)
iasva.sv <- iasva.res$sv
plot(iasva.sv[, 1], iasva.sv[, 2], xlab = "SV1", ylab = "SV2")
cell_type <- as.factor(iasva.sv[, 1] > -0.1) 
levels(cell_type) <- c("Cell1", "Cell2")
table(cell_type)

# We identified 6 outlier cells based on SV1 that are marked in red
pairs(iasva.sv, main = "IA-SVA", pch = 21, col = color.vec[cell_type],
      bg = color.vec[cell_type], oma = c(4,4,6,12))
legend("right", levels(cell_type), fill = color.vec, bty = "n")
plot(iasva.sv[, 1:2], main = "IA-SVA", pch = 21, xlab = "SV1", ylab = "SV2",
     col = color.vec[cell_type], bg = color.vec[cell_type])
cor(geo_lib_size, iasva.sv[, 1])
corrplot(cor(iasva.sv))
```

As shown in the above figure, SV1 clearly separates alpha cells into two
groups: 6 outlier cells (marked in red) and the rest of the alpha cells
(marked in blue). 
## Find marker genes for the detected heterogeneity (SV1).
Here, using the find_markers() function we find marker genes that are
significantly associated with SV1 (multiple testing adjusted p-value < 0.05,
default significance cutoff, and R-squared value > 0.3, default R-squared cutoff).  
```{r find_markers, echo=TRUE, fig.width=7, fig.height=14}
marker.counts <- find_markers(summ_exp, as.matrix(iasva.sv[,1]))
nrow(marker.counts)
rownames(marker.counts)
anno.col <- data.frame(cell_type = cell_type)
rownames(anno.col) <- colnames(marker.counts)
head(anno.col)
pheatmap(log(marker.counts + 1), show_colnames = FALSE,
         clustering_method = "ward.D2", cutree_cols = 2,
         annotation_col = anno.col)
```

## Run tSNE to detect the hidden heterogeneity.
For comparison purposes, we applied tSNE on read counts of all genes to
identify the hidden heterogeneity. We used the Rtsne R package with default settings.

```{r run_tsne, echo=TRUE, fig.width=7, fig.height=7}
set.seed(100)
tsne.res <- Rtsne(t(lcounts), dims = 2)

plot(tsne.res$Y, main = "tSNE", xlab = "tSNE Dim1", ylab = "tSNE Dim2",
     pch = 21, col = color.vec[cell_type], bg = color.vec[cell_type],
     oma = c(4, 4, 6, 12))
legend("bottomright", levels(cell_type), fill = color.vec, bty = "n")
```

As shown above, tSNE fails to detect the outlier cells that are identified by
IA-SVA when all genes are used. Same color coding is used as above. 

## Run tSNE post IA-SVA analyses, i.e., run tSNE on marker genes associated with SV1 as detected by IA-SVA. 
Here, we apply tSNE on the marker genes for SV1 obtained from IA-SVA
```{r run_tsne_post_iasva, echo=TRUE, fig.width=7, fig.height=7}
set.seed(100)
tsne.res <- Rtsne(unique(t(log(marker.counts + 1))), dims = 2)
plot(tsne.res$Y, main = "tSNE post IA-SVA", xlab = "tSNE Dim1",
     ylab = "tSNE Dim2", pch = 21, col = color.vec[cell_type],
     bg = color.vec[cell_type], oma = c(4, 4, 6, 12))
legend("bottomright", levels(cell_type), fill = color.vec, bty = "n")
```

tSNE using SV1 marker genes better seperate these ourlier cells. 
This analyses suggest that gene selection using IA-SVA combined
with tSNE analyses can be a powerful way to detect rare cells introducing
variability in the single cell gene expression data. 

## Using a faster implementation of IA-SVA (fast_iasva)
Here, we run a faster implementation of IA-SVA using the same known
factors (patient_id and geo_lib_size) as demonstrated above. This function
is useful when working with particularly large datasets.
```{r fast_iasva, echo=TRUE}
iasva.res <- fast_iasva(summ_exp, mod[, -1], num.sv = 5)
```

## Tuning parameters for IA-SVA
The R-squared thresholds used to identify marker genes (find_markers)
can greatly influence the 1) number of marker genes identified and 2)
the quality of clustering results. With the study_R2() function,
users can visualize how different R-squared thresholds influence both factors.

```{r study_r2, echo=TRUE, fig.width=7, fig.height=7}
study_res <- study_R2(summ_exp, iasva.sv)
```
This function produces a plot of the number of genes selected vs. the 
cluster quality (average silhouette score) for different R-squared values.

## Session Info
```{r session_info, echo=TRUE}
sessionInfo()
```
