paul-dim-red.Rmd

---
title: "Paul Gene Expression Data: Dimension Reduction"
author: "Michelle Li and Gianna Wu with Professor de Pillis and Dr. Park"
date: "7/20/2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(destiny)
library(Biobase)
library(tsne)
library(rgl)
options(rgl.useNULL = TRUE)
r3dDefaults$windowRect <- c(0,0,1000,1000)
```

```{r load-paul, include=FALSE}
## load .Rdata file into environment so you have all calculations done (located in our Drive)
## if you have the .RData file to load, use eval=FALSE for pca, diffmap, tsne
load("~/GitHub/cell-diff-reu-2018/paul-rmd-data.RData")
```

----

Analyzing Paul gene expression data using the following dimension reduction techniques in R: PCA (handwritten version), diffusion mapping (destiny package), and t-SNE (tsne package). The original Paul data comes in as a matrix of 8716 genes x 2730 cells.  

Paul examines stem cell differentiation in the bone marrow. Using single-cell RNA sequencing they examine the heterogeneity of myeloid progenitor cells to better understand and define the differentiation of hematopoietic progenitors into distinct
lineages. They identify 19 transcriptionally homogeneous subpopulations of cells, which we use for our coloring of plots.  

```{r prepare-paul, eval=FALSE, include=FALSE}
## read in data 
load("~/GitHub/cell-diff-reu-2018/data/Paul_Cell_MARSseq_GSE72857.RData")

## create factor for grouping
groupsf <- factor(cluster.id) # paul data
nclusters <- nlevels(groupsf) # number of groups to cluster
colors <- topo.colors(nclusters) # choose colors
```

```{r heatmap-grover, eval=FALSE, include=FALSE}
hmap <- heatmap(data, Rowv=NA, Colv=NA)
```

## PCA (Handwritten)

```{r pca-paul, eval=FALSE, include=FALSE}
## PCA DIMENSIONS:
## cells = rows, genes = columns 
## Paul: use TRANSPOSE of data matrix (want cells x genes)
data <- t(data)
M1centered <- scale(data, center=TRUE, scale=FALSE) # subtracts column means from each column

## svd
svd <- svd(M1centered, nu=3, nv=3) # only first 3 left/right singular vectors for efficiency
U <- svd$u # dim: nrowx3
V <- svd$v # dim: ncolx3
# Note: U = M1centered %*% V %*% S^-1

## calculate scores
# weights <- t(M1centered) %*% U # dim: ncolx3
scores <- M1centered %*% V # dim: nrowx3
```

#### Plots of scores 
Scores were calculated by: scores <- data %*% V. V are the right singular vectors from SVD whose columns contain the eigenvectors of the covariance matrix. Scores give the coordinates of the individuals on the principal components (coefficients of the linear combination of original variables for each principal component). 

2D plot of scores   
```{r pca-plot2dscores}
plot(x=scores[,1], y=scores[,2], 
     xlab="PC1 scores", ylab="PC2 scores", 
     col = colors[groupsf], 
     pch=20, 
     main = "PCA Scores, Paul")
legend("bottomright", inset = c(0,0), 
       legend = levels(groupsf), 
       pch = 20, 
       col = colors, 
       ncol=3, cex=0.4, pt.cex=1)
```

3D plot of scores 
```{r pca-plot3dscores}
plot3d(x=scores[,1], y=scores[,2], z=scores[,3], 
       xlab="PC1 scores", ylab="PC2 scores", zlab="PC3 scores", 
       col = colors[groupsf], 
       pch=20, 
       main = "PCA Scores, Paul")
legend3d("bottomright", inset=c(0.1,0.1),
         legend = levels(groupsf), 
         pch = 20, 
         col = colors, 
         ncol=3, cex=1.5)
rglwidget()
```

#### Plots of principal components

2D plot of PCs (eigenvectors of covariance matrix)  
```{r pca-plot2dpcs}
## plot PCs 2d (eigenvectors of covariance matrix)
plot(x=V[,1], y=V[,2], 
     xlab="PC1", ylab="PC2", 
     col = colors[groupsf], 
     pch=20, 
     main = "Principal Components, Paul")
legend("bottomright", inset = c(0,0),
       legend = levels(groupsf),
       pch = 20,
       col = colors,
       ncol=3, cex=0.4, pt.cex=1)
```

3D plot of PCs (eigenvectors of covariance matrix)  
```{r pca-plot3dpcs}
## plot PCs 3d
plot3d(x=V[,1], y=V[,2], z=V[,3], 
       xlab="PC1", ylab="PC2", zlab="PC3", 
       col = colors[groupsf], 
       pch=20, 
       main = "Principal Components, Paul")
legend3d("bottomright", inset=c(0.1,0.1),
         legend = levels(groupsf),
         pch = 20, 
         col = colors, 
         ncol=3, cex=1.5)
rglwidget()
```

----

## Diffusion Mapping (destiny)

Diffusion mapping using destiny package.  
```{r diffmap-paul, eval=FALSE, include=FALSE}
## convert to expression set
## EXPRESSION SET DIMENSIONS: cells = columns ; genes = rows
## Paul needs the original formatting (we just took a transpose, must transpose back)
data <- t(data)
dataES <- ExpressionSet(assayData=data) # you must include the 'assayData=' for this to come out correctly

## diffusion map
dm <- DiffusionMap(dataES)

## sigmas
sigmas <- find_sigmas(dataES) 
```

2D Diffusion Map  
```{r diffmap-plot2d}
## plot 2d
plot(x=eigenvectors(dm)[,1], y=eigenvectors(dm)[,2],
     xlab="DC1", ylab="DC2",
     col = colors[groupsf],
     pch=20, 
     main = "Diffusion Map, Paul")
legend("bottomright", inset=c(0,0), 
       legend=levels(groupsf), 
       col=colors, 
       pch=20, 
       ncol=3, cex=0.4, pt.cex=1) 
```

3D Diffusion Map  
```{r diffmap-plot3d}
## plot 3d
plot3d(x=eigenvectors(dm)[,1], y=eigenvectors(dm)[,2], z=eigenvectors(dm)[,3],
       xlab="DC1", ylab="DC2", zlab="DC3",
       col = colors[groupsf],
       pch=20, 
       main = "Diffusion Map, Paul")
legend3d("bottomright", inset=c(0.1,0.1),
         legend = levels(groupsf),
         pch = 20, 
         col = colors, 
         ncol=3, cex=1.5)
rglwidget()
```

Sigmas  
```{r diffmap-sigmas}
plot(sigmas)
```

----

## t-SNE (tsne package)

t-SNE using tsne package (from jdonaldson, linked on van der Maaten's website).

```{r tsne-paul, eval=FALSE, include=FALSE}
## already read in data and converted to matrix
# Note: can take very, very long with this dataset
# t-SNE DIMENSIONS: rows=cells, columns=genes
# Paul data needs transpose to have correct dimensions
data <- t(data)
ydata3d <- tsne(data, k=3, max_iter = 200) # ydata = matrix(rnorm(k * n),n) 
```

2D t-SNE plot  
```{r tsne-plot2d}
plot(x=ydata3d[,1], y=ydata3d[,2], 
     xlab="", ylab="",
     col = colors[groupsf], 
     pch=20,
     main = "t-SNE, Paul") 
legend("bottomleft", inset = c(0,0), 
       legend = levels(groupsf), 
       pch = 20, 
       col = colors, 
       ncol=3, cex=0.4, pt.cex=1) 
```

3D t-SNE plot  
```{r tsne-plot3d}
plot3d(x=ydata3d[,1], y=ydata3d[,2], z=ydata3d[,3], 
       xlab="", ylab="", zlab="",
       col=colors[groupsf], 
       pch=20,
       main = "t-SNE, Paul")
legend3d("bottomright", inset=c(0.1,0.1),
         legend = levels(groupsf),  
         pch = 20, 
         col = colors, 
         ncol=3, cex=1.5)
rglwidget()
```