UCell_demo.Rmd

---
title: UCell demo
author: 
- Massimo Andreatta^[massimo.andreatta@unil.ch]
- Santiago Carmona^[santiago.carmona@unil.ch]
date: "15/02/2021"
knit: (function(input_file, encoding) {
  out_dir <- 'docs';
  rmarkdown::render(input_file,
 encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'index.html'))})
#output: html_notebook
---


```{r message=F, warning=F}
#renv::restore()

BiocManager::install("UCell")
#remotes::install_github("carmonalab/UCell")  #development version

library(UCell)
```

#Object at: https://drive.switch.ch/index.php/s/3kM5PQ0tQaG6d6A
```{r}
library(Seurat)
pbmc.Tcell <- readRDS("data/pbmc_multimodal.downsampled20k.Tcell.seurat.RNA.rds")
DimPlot(object = pbmc.Tcell, reduction = "wnn.umap", group.by = "celltype.l2", label = TRUE, label.size = 3, repel = TRUE)
```
# Basic functionalities
```{r}
data.matrix <- pbmc.Tcell@assays$RNA@data

dim(data.matrix)
data.matrix[1:5,1:5]
```

Calculate on the fly
```{r}
set.seed(123)
#basic.sign <- list( Tcell_signature = c("CD2","CD3E","CD3D"), Myeloid_signature = c("SPI1","FCER1G","CSF1R"))
basic.sign <- list( Tcell_signature = c("CD2","CD3E","CD3D"), Myeloid_signature = c("SPI1","FCER1G","CSF1R","LCK-"))


scores <- ScoreSignatures_UCell(data.matrix,features=basic.sign)
scores[1:5,]
```

Show distribution of predicted scores
```{r}
library(reshape2)
library(ggplot2)
melted <- melt(scores)
colnames(melted) <- c("Cell","Signature","Uscore")
p <- ggplot(melted, aes(x=Signature, y=Uscore)) + geom_violin(aes(fill=Signature), scale = "width") +
       geom_boxplot(width=0.1) + theme_bw()
p
```

Pre-compute ranks
```{r}
set.seed(123)
ranks <- StoreRankings_UCell(data.matrix)
ranks[1:5,1:5]
```

Now quickly evaluate signatures from pre-stored ranks
```{r}
set.seed(123)
scores.2 <- ScoreSignatures_UCell(data.matrix, features=basic.sign, precalc.ranks = ranks)

scores.2[1:5,]

plot(scores[,1], scores.2[,1])


new.sign <- list( CD4_Tcell = c("CD4","CD40LG"), CD8_Tcell = c("CD8A","CD8B"), Treg = c("FOXP3","IL2RA"))
scores.new <- ScoreSignatures_UCell(data.matrix, features=new.sign, precalc.ranks = ranks)

scores.new[1:5,]
melted <- melt(scores.new)
colnames(melted) <- c("Cell","Signature","Uscore")
p <- ggplot(melted, aes(x=Signature, y=Uscore)) + geom_violin(aes(fill=Signature), scale = "width") +
       geom_boxplot(width=0.1) + theme_bw()
p
```


Run in parallel with multiple cores
```{r}
scores <- ScoreSignatures_UCell(data.matrix, features=basic.sign, ncores=4)

melted <- melt(scores)
colnames(melted) <- c("Cell","Signature","Uscore")
p <- ggplot(melted, aes(x=Signature, y=Uscore)) + geom_violin(aes(fill=Signature), scale = "width") +
       geom_boxplot(width=0.1) + theme_bw()
p
```


# Easy interaction with Seurat

Define human T cell level 2 signatues
```{r, eval=T}
signatures.T <- list()
signatures.T$Tcell_CD4 <- c("CD4","CD40LG")
signatures.T$Tcell_CD8 <- c("CD8A","CD8B")
signatures.T$Tcell_Treg <- c("FOXP3","IL2RA")
signatures.T$Tcell_MAIT <- c("KLRB1", "SLC4A10", "NCR3")
signatures.T$Tcell_gd <- c("TRDC", "TRGC1", "TRGC2", "TRDV1","TRAC-","TRBC1-","TRBC2-")
signatures.T$Tcell_NK <- c("FGFBP2", "SPON2", "KLRF1", "FCGR3A", "KLRD1", "TRDC","CD3E-","CD3G-")
```



```{r}
Idents(pbmc.Tcell) <- "celltype.l1"
pbmc.Tcell.CD8 <- subset(pbmc.Tcell, idents = c("CD8 T"))
table(pbmc.Tcell.CD8@active.ident)
DimPlot(object = pbmc.Tcell.CD8, reduction = "wnn.umap", group.by = "celltype.l2", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
```

AddModuleScore from Seurat is very fast, but the score is highly dependent on the composition of the dataset
```{r}
gene.sets <- list(Tcell_CD8=c("CD8A","CD8B"))

pbmc.Tcell.CD8 <- AddModuleScore(pbmc.Tcell.CD8, features = gene.sets, name="Tcell_CD8_Seurat", seed=123)
p1 <- VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_Seurat1")

pbmc.Tcell <- AddModuleScore(pbmc.Tcell, features = gene.sets, name="Tcell_CD8_Seurat", seed=123)
p2 <- VlnPlot(pbmc.Tcell, features = "Tcell_CD8_Seurat1")

p1 | p2

```

UCell score is based on gene rankings and therefore is not affected by the composition of the query dataset
```{r}
pbmc.Tcell.CD8 <- AddModuleScore_UCell(pbmc.Tcell.CD8, features = gene.sets)
p1 <- VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_UCell")

pbmc.Tcell <- AddModuleScore_UCell(pbmc.Tcell, features = gene.sets)
p2 <- VlnPlot(pbmc.Tcell, features = "Tcell_CD8_UCell")

p1 | p2
```

Evaluate more signatures for all flavors of T cells
```{r fig.height=3}
pbmc.Tcell.ss  <- AddModuleScore_UCell(pbmc.Tcell, features = signatures.T)
sign.names <- paste0(names(signatures.T),"_UCell")

VlnPlot(pbmc.Tcell.ss, features = sign.names)
VlnPlot(pbmc.Tcell.ss, features = sign.names, group.by = "celltype.l2")
```

How do signatures compare to original annotations
```{r fig.height=3}
Idents(pbmc.Tcell.ss) <- "celltype.l2"
DimPlot(object = pbmc.Tcell.ss, reduction = "wnn.umap", group.by = "celltype.l2", label.size = 3, repel = TRUE, label = T)
FeaturePlot(pbmc.Tcell.ss, reduction = "wnn.umap", features = sign.names, ncol=3, order=T)
```

See behavior with missing or non-existing genes
```{r}
pbmc.Tcell.CD8 <- AddModuleScore_UCell(pbmc.Tcell.CD8, features = list(Tcell_CD8=c("CD8A","CD8B","notagene")))
p1 <- VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_UCell")

pbmc.Tcell.CD8 <- AddModuleScore_UCell(pbmc.Tcell.CD8, features = list(Tcell_CD8=c("CD8A","CD8B","notagene","notagene2","notagene3")))
p2 <- VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_UCell")

pbmc.Tcell.CD8 <- AddModuleScore_UCell(pbmc.Tcell.CD8, features = list(Tcell_CD8=c("notagene","notagene2","notagene3")))
p3 <- VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_UCell")

p1 | p2 | p3

```


Recalculate with the Cd8 T cell simple signature
```{r}
gene.sets <- list(Tcell_CD8=c("CD8A","CD8B"))
pbmc.Tcell.CD8 <- AddModuleScore_UCell(pbmc.Tcell.CD8, features = gene.sets)
VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_UCell")


plot(pbmc.Tcell.CD8$Tcell_CD8_Seurat1, pbmc.Tcell.CD8$Tcell_CD8_UCell)
plot(pbmc.Tcell$Tcell_CD8_Seurat1, pbmc.Tcell$Tcell_CD8_UCell)

```


Compare to AUCell
```{r}
library(AUCell)
cells_rankings <- AUCell_buildRankings(pbmc.Tcell@assays$RNA@data,nCores = 1, plotStats = F)
cells_AUC <- AUCell_calcAUC(signatures.T, cells_rankings, aucMaxRank=1000)
cells_AUC.num <- as.data.frame(t(getAUC(cells_AUC)))
colnames(cells_AUC.num) <- paste0(colnames(cells_AUC.num), "_AUCell")

pbmc.Tcell <- AddMetaData(pbmc.Tcell, cells_AUC.num)


cells_rankings.cd8 <- AUCell_buildRankings(pbmc.Tcell.CD8@assays$RNA@data,nCores = 1, plotStats = T)
cells_AUC.cd8 <- AUCell_calcAUC(signatures.T, cells_rankings.cd8, aucMaxRank=1000)
cells_AUC.num.cd8 <- as.data.frame(t(getAUC(cells_AUC.cd8)))
colnames(cells_AUC.num.cd8) <- paste0(colnames(cells_AUC.num.cd8), "_AUCell")

pbmc.Tcell.CD8 <- AddMetaData(pbmc.Tcell.CD8, cells_AUC.num.cd8)


p1 <- VlnPlot(pbmc.Tcell, features = "Tcell_CD8_AUCell")
p2 <- VlnPlot(pbmc.Tcell.CD8, features = "Tcell_CD8_AUCell")
p1 | p2
```



```{r}
#Seurat ModuleScore vs UCell
plot(pbmc.Tcell.CD8$Tcell_CD8_Seurat1, pbmc.Tcell.CD8$Tcell_CD8_UCell)

#AUCell vs UCell
plot(pbmc.Tcell.CD8$Tcell_CD8_AUCell, pbmc.Tcell.CD8$Tcell_CD8_UCell)
```



## Test ranks storing functionality

```{r eval=F}
t1 <- Sys.time()
pbmc.Tcell.UcellPrecomp <- AddModuleScore_UCell(pbmc.Tcell, features = gene.sets, storeRanks = T)
t2 <- Sys.time()
t2 - t1

VlnPlot(pbmc.Tcell.UcellPrecomp, features = "Tcell_CD8_UCell")
```

Check that scores calculate with rank pre-calculation are the same as those calculated "on the fly"
```{r eval=F}
t1 <- Sys.time()
pbmc.Tcell.UcellPrecomp <- AddModuleScore_UCell(pbmc.Tcell.UcellPrecomp, features = list(Tcell_CD8_copy=c("CD8A","CD8B")))
t2 <- Sys.time()
t2 - t1

VlnPlot(pbmc.Tcell.UcellPrecomp, features = "Tcell_CD8_copy_UCell")

plot(pbmc.Tcell.UcellPrecomp$Tcell_CD8_UCell,pbmc.Tcell.UcellPrecomp$Tcell_CD8_copy_UCell)

t1 <- Sys.time()
pbmc.Tcell.UcellPrecomp <- AddModuleScore_UCell(pbmc.Tcell.UcellPrecomp, features = signatures.T)
t2 <- Sys.time()
t2 - t1

VlnPlot(pbmc.Tcell.UcellPrecomp, features = paste0(names(signatures.T),"_UCell"))

```

Gene ranks are stored here:
```{r eval=F}
head(pbmc.Tcell.UcellPrecomp@assays$UCellRanks@data[,1:10])
head(pbmc.Tcell.UcellPrecomp@assays$UCellRanks@data["CD8A",1:10])
hist(as.numeric(pbmc.Tcell.UcellPrecomp@assays$UCellRanks@data["CD8A",]))
```


# Test escape enrichIT
http://www.bioconductor.org/packages/release/bioc/vignettes/escape/inst/doc/vignette.html#3_Enrichment

```{r}
devtools::install_github("ncborcherding/escape")

library(escape)
#library(dittoSeq)
#library(SingleCellExperiment)
#pbmc.sce <- as.SingleCellExperiment(pbmc.Tcell)
```

```{r}
library(GSEABase)
#check format for gene sets
GS <- getGeneSets(library = "H")

gsets <- list()
for (sig in names(signatures.T)) {
   gsets[[sig]] <- GeneSet(signatures.T[[sig]], setName=sig)
}
gsc <- GeneSetCollection(gsets)

esc.res <- enrichIt(pbmc.Tcell, gene.sets = gsc, groups = 1000)

```

Very slow, and gives scores in strange multimodal distributions...
```{r}
colnames(esc.res) <- paste0(colnames(esc.res), "_enrichIT")

pbmc.Tcell <- AddMetaData(pbmc.Tcell, esc.res)
VlnPlot(pbmc.Tcell, features = colnames(esc.res))
```

Export small subset as sample data for the package
```{r}
Idents(pbmc.Tcell) <- "celltype.l1"
pbmc.ds <- subset(pbmc.Tcell, cells=WhichCells(pbmc.Tcell,downsample=200))
sample.matrix <- pbmc.ds@assays$RNA@data
save(sample.matrix, file="sample.matrix.RData")
```


# Test ties methods

```{r}
set.seed(123)
basic.sign <- list( Tcell_signature = c("CD2","CD3E","CD3D"), Myeloid_signature = c("SPI1","FCER1G","CSF1R"))

system.time(scores <- ScoreSignatures_UCell(data.matrix,features=basic.sign, ties.method = "random"))
scores[1:5,]

#Using average for rank ties
system.time(scores <- ScoreSignatures_UCell(data.matrix,features=basic.sign, ties.method = "average"))
scores[1:5,]

scores.rand <- list()
scores.ave <- list()
for (seed in c(1:3)) {
   set.seed(seed)
   t1 <- system.time(scores.rand[[seed]] <- ScoreSignatures_UCell(data.matrix,features=basic.sign, ties.method = "random"))
   t2 <- system.time(scores.ave[[seed]] <- ScoreSignatures_UCell(data.matrix,features=basic.sign, ties.method = "ave"))
   print(t1)
   print(t2)
}

##Reproducibility
#Rand vs rand
plot(scores.rand[[1]][,1], scores.rand[[2]][,1])
plot(scores.rand[[1]][,1], scores.rand[[3]][,1])

#ave vs ave
plot(scores.ave[[1]][,1], scores.ave[[2]][,1])
plot(scores.ave[[1]][,1], scores.ave[[3]][,1])

#Rand vs. ave
plot(scores.rand[[1]][,1], scores.ave[[1]][,1])
plot(scores.rand[[2]][,1], scores.ave[[2]][,1])
```