clusterW.Rmd

---
title: "Cluster cells using seurat on W - 10x genomics"
author: "Fanny Perraudeau"
date: "`r Sys.Date()`"
output: 
html_document: 
fig_height: 7
fig_width: 7
toc: yes
code_folding: hide
toc_float: yes
---

```{r options, echo=FALSE, results="hide",mesasge=FALSE, error=FALSE, include=FALSE, autodep=TRUE}
knitr::opts_chunk$set(fig.align="center", cache=TRUE, error=FALSE, message=FALSE, warning=TRUE)
library(zinbwave)
library(Seurat)
library(Rtsne)
set.seed(63738)
```

We want to perform clustering on dimensionality reduced matrix W estimated using zinbwave. We use Seurat method to perform clustering but tweak the function to be able to into the W matrix instead of the reduced dimension matrix from PCA.

# Data
```{r}
# SingleCellExperiment object
load("../data/core.rda")
core = core[rowSums(assay(core)) > 0, ]
colData(core)$seurat = factor(colData(core)$seurat)
core = core[rowData(core)$seuratVarGenes, ]
dim(core)
```

```{r}
# zinb 
# tried K=10,25 - same as PCA in seurat tuto
load('zinb_k20.rda')
W = getW(zinb)
dim(W)
```

```{r}
# seurat object
load('../data/pbmc3k.rda')
pbmc
seurat_clus <- pbmc@ident
names(seurat_clus) <- paste0(names(pbmc@ident), "-1")
pbmc_seurat = pbmc
current.cluster.ids <- c(0, 1, 2, 3, 4, 5, 6, 7, 8)
new.cluster.ids <- c("CD4 T cells memory","CD4 T cells naive", "CD14+ Monocytes", "B cells", "CD8 T cells", "FCGR3A+ Monocytes", "NK cells", "Dendritic cells", "Megakaryocytes")
pbmc_seurat@ident <- plyr::mapvalues(x = pbmc_seurat@ident, from = current.cluster.ids, to = new.cluster.ids)
```


# Seurat clustering on W

Use Davide Risso's code to try Seurat clustering on W. We need to change the code in Seurat so that W matrix from zinbwave is used instead of the pca coordinates.

```{r seuratClus}
library(FNN)
library(dplyr)
library(Matrix)

## Build SNN
## Keep default parameters from seurat R function FindClusters()
k.param = 30
k.scale = 25
prune.SNN = 1/15
modularity.fxn = 1
resolution = 1.8
n.start = 100
n.iter = 10
algorithm = 1
random.seed = 0
n.cells = nSamples(zinb)
data.use <- W

# KNN
my.knn <- FNN::get.knn(as.matrix(data.use),
                       k = min(k.scale * k.param, n.cells - 1))
nn.ranked <- cbind(1:n.cells, my.knn$nn.index[, 1:(k.param-1)])
nn.large <- my.knn$nn.index
# SNN
w <- Seurat:::CalcSNNSparse(colnames(pbmc@data), k.param,
                            nn.large, nn.ranked, 
                            prune.SNN = prune.SNN, 
                            print.output = FALSE)
## Run modularity clustering
SNN.use <- w
pbmc <- Seurat:::RunModularityClustering(pbmc, SNN.use, 
                                         modularity = modularity.fxn,
                                         r = resolution,
                                         algorithm = algorithm,
                                         n.start = n.start, 
                                         n.iter = n.iter, 
                                         random.seed = random.seed, 
                                         print.output = FALSE,
                                         temp.file.location = NULL)
pbmc <- Seurat:::GroupSingletons(pbmc, SNN.use)
name <- paste("res.", resolution, sep = "")
pbmc <- StashIdent(pbmc, name)
#```
#
#Compare Seurat and ours clusterings.
#```{r}
our_clus <- pbmc@ident
names(our_clus) <- paste0(names(pbmc@ident), "-1")
table(seurat_clus, our_clus)
```
```{r}
save(pbmc, file = 'pmbc.rda')
save(pbmc_seurat, file = 'pmbc_seurat.rda')
```

```{r}
set.seed(9374)
zinb_tsne <- Rtsne(data.use, pca = FALSE, perplexity = 200)
zinb_tsne <- as.data.frame(zinb_tsne$Y)
colnames(zinb_tsne) = c('tSNE_1', 'tSNE_2')
```

```{r}
library(plyr)
current <- c(0, 1, 2, 3, 4, 5, 6, 7, 8)
new <- c("CD4 T cells memory", "CD4 T cells naive", "CD14+ Monocytes", "B cells", "CD8 T cells", "FCGR3A+ Monocytes", "NK cells", "Dendritic cells", "Megakaryocytes")
new_seurat <- plyr::mapvalues(x = seurat_clus,
                              from = current, to = new)
```

```{r}
library(ggplot2)
seuratTheme = theme_bw() + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

p1 <- ggplot(zinb_tsne, aes(x = tSNE_1, y = tSNE_2)) + 
  geom_point(mapping = aes(colour = factor(x = our_clus))) +
  scale_color_brewer(palette="Set1") + seuratTheme 

p2 <- ggplot(zinb_tsne, aes(x = tSNE_1, y = tSNE_2)) + 
  geom_point(mapping = aes(colour = factor(x = new_seurat))) +
  guides(col=guide_legend(title="PCA")) + seuratTheme 

p3 <- TSNEPlot(object = pbmc, do.return = TRUE, no.legend = FALSE,
               do.label = FALSE)
p4 <- TSNEPlot(object = pbmc_seurat, do.return = TRUE, no.legend = FALSE, do.label = FALSE, pt.size = 0.5)

legend_seurat <- get_legend(p2 + guides(col=guide_legend(title="PCA")))
legend_ours <- get_legend(p1 + guides(col=guide_legend(title="W")))
legendall <- plot_grid(legend_ours, legend_seurat, ncol = 1)

pall = plot_grid(p1 + theme(legend.position = "none") ,
                 p2 + theme(legend.position = "none"),
                 p3 + scale_color_brewer(palette="Set1") +
                   theme(legend.position = "none"),
                 p4 + theme(legend.position = "none"))
plot_grid(pall, legendall, rel_widths = c(3, 1))
```

```{r tenxcaseW,fig.height=3,fig.width=10}
library(cowplot)
plot_grid(p3 + theme(legend.position = "none") +
            scale_color_brewer(palette="Set1"),
          p1 + theme(legend.position = "none"),
          p2 + theme(legend.position = "none"),
          labels = c("a", "b", "c"), ncol = 3, nrow = 1,
          align = "h")
```

```{r}
size=20
```

```{r tenxcaseWW}
p1 <- ggplot(zinb_tsne, aes(x = tSNE_1, y = tSNE_2)) + 
  geom_point(mapping = aes(colour = factor(x = our_clus))) +
  scale_color_brewer(palette="Set1") + seuratTheme 
p1 + labs(color='Clustering with W') +
            theme(text = element_text(size = size)) +
            guides(colour = guide_legend(override.aes = list(size = 5)))
```

```{r tenxcaseWPCA}
p2 <- ggplot(zinb_tsne, aes(x = tSNE_1, y = tSNE_2)) + 
  geom_point(mapping = aes(colour = factor(x = new_seurat))) +
  guides(col=guide_legend(title="PCA")) + seuratTheme 
p2 + labs(color='Clustering with PCA') +
  guides(colour = guide_legend(override.aes = list(size = 5)))+
  theme(text = element_text(size=size))
```

```{r tenxcasePCAW}
p3 <- TSNEPlot(object = pbmc, do.return = TRUE, no.legend = FALSE, do.label = FALSE)
p3 + scale_color_brewer(palette="Set1") +
  guides(colour = guide_legend(override.aes = list(size = 5))) +
  theme(text = element_text(size = size), legend.title = element_text()) +
  labs(colour = 'Clustering with W')
```

```{r tenxcasePCAPCA}
p4 <- TSNEPlot(object = pbmc_seurat, do.return = TRUE, no.legend = FALSE,
               do.label = FALSE)
p4 + labs(color='Clustering with PCA') +
  guides(colour = guide_legend(override.aes = list(size = 5))) +
  theme(text = element_text(size=size), legend.title = element_text()) +
  labs(colour = 'Clustering with PCA')
```

# Seurat DE on 2 clusters for CD4+ cells

```{r}
tcell.markers <- FindMarkers(object = pbmc,
                             ident.1 = 0,
                             ident.2 = 5,
                             logfc.threshold = 0,
                             genes.use = rownames(pbmc@raw.data),
                             thresh.use = 0,
                             test.use = "wilcox",
                             min.pct = 0, 
                             min.diff.pct = -Inf,
                             print.bar = TRUE,
                             only.pos = FALSE,
                             max.cells.per.ident = Inf, 
                             random.seed = 1,
                             latent.vars = "nUMI",
                             min.cells = 0 )
```


```{r restore}
library(dplyr)
head(tcell.markers, 2)
length(unique(rownames(pbmc@raw.data)))
length(unique(rownames(tcell.markers)))
tcell.markers$gene = rownames(tcell.markers)
topgenes = tcell.markers %>% top_n(20, -p_val)
topgenes
```

```{r}
FeaturePlot(object = pbmc,
            features.plot = c("NKG7", "S100A4"), 
            cols.use = c("green", "blue"))
```

```{r}
write.csv(tcell.markers, file = 'tcellmarkers_seurat.csv')
```

# Save new SingleCellExeriment object
```{r}
colData(core)$ours = our_clus
save(core, file = 'coreAfterClustering.rda')
```

# sessionInfo

```{r}
sessionInfo()
```