---
author: "Åsa Björklund  &  Paulo Czarnewski"
date: '`r format(Sys.Date(), "%B %d, %Y")`'
output:
  html_document:
    self_contained: true
    highlight: tango
    df_print: paged
    toc: yes
    toc_float:
      collapsed: false
      smooth_scroll: true
    toc_depth: 3
    keep_md: yes
    fig_caption: true
  html_notebook:
    self_contained: true
    highlight: tango
    df_print: paged
    toc: yes
    toc_float:
      collapsed: false
      smooth_scroll: true
    toc_depth: 3
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(message=FALSE, warning=FALSE, result='hold',fig.width=12,tidy=TRUE)
knitr::opts_knit$set(progress=TRUE,verbose=TRUE)
```
<style>
h1, .h1, h2, .h2, h3, .h3, h4, .h4 { margin-top: 50px }
p.caption {font-size: 0.9em;font-style: italic;color: grey;margin-right: 10%;margin-left: 10%;text-align: justify}
</style>

## Celltype prediction
***

 Celltype prediction can either be performed on indiviudal cells where each cell gets a predicted celltype label, or on the level of clusters. All methods are based on similarity to other datasets, single cell or sorted bulk RNAseq, or uses known marker genes for each celltype.

We will select one sample from the Covid data, `ctrl_13` and predict celltype by cell on that sample.

Some methods will predict a celltype to each cell based on what it is most similar to even if the celltype of that cell is not included in the reference. Other methods include an uncertainty so that cells with low similarity scores will be unclassified.
There are multiple different methods to predict celltypes, here we will just cover a few of those. 

Here we will use a reference PBMC dataset from the `scPred` package which is already a Seurat object with counts. And we will test classification based on label transfer using the function `TransferData` in the Seurat package and the `scPred` method. Finally we will use gene set enrichment predict celltype based on the DEGs of each cluster. 

# Load and process data
## Covid-19 data
First, lets load required libraries and the saved object from the clustering step. Subset for one patient.

```{r}
suppressPackageStartupMessages({
  library(Seurat)
  library(dplyr)
  library(cowplot)
  library(ggplot2)
  library(pheatmap)
  library(rafalib)
  library(scPred)
})
```

```{r}
#load the data and select 'ctrl_13` sample
alldata <- readRDS("data/results/covid_qc_dr_int_cl.rds")
ctrl = alldata[, alldata$orig.ident == 'ctrl_13']

# set active assay to RNA and remove the CCA assay
ctrl@active.assay = 'RNA' 
ctrl[['CCA']] = NULL
ctrl
```


## Reference data
Then, load the reference dataset with annotated labels. Also, run all steps of the normal analysis pipeline with normalizaiton, variable gene selection, scaling and dimensionality reduction.

```{r}
reference <- scPred::pbmc_1

reference
```


## Rerun analysis pipeline
Here, we will run all the steps that we did in previous labs in one go using the `magittr` package with the pipe-operator `%>%`.

```{r}
reference <- reference %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = F) %>%
  RunUMAP(dims = 1:30)
```

```{r, fig.width=5}
DimPlot(reference, group.by = "cell_type", label = TRUE, repel = TRUE) + NoAxes()
```


Run all steps of the analysis for the ctrl sample as well. Use the clustering from the integration lab with resolution 0.5.

```{r}
#Set the identity as louvain with resolution 0.3
ctrl <- SetIdent(ctrl, value = "CCA_snn_res.0.5")
  
ctrl <- ctrl %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = F) %>%
  RunUMAP(dims = 1:30)
```

```{r, fig.width=5}
DimPlot(ctrl,  label = TRUE, repel = TRUE) + NoAxes()
```


# Seurat label transfer
First we will run label transfer using a similar method as in the integration exercise. But, instad of CCA the default for the 'FindTransferAnchors` function is to use "pcaproject", e.g. the query dataset is projected onto the PCA of the reference dataset. Then, the labels of the reference data are predicted.

```{r}
transfer.anchors <- FindTransferAnchors(reference = reference, query = ctrl, 
    dims = 1:30)
predictions <- TransferData(anchorset = transfer.anchors, refdata = reference$cell_type, 
    dims = 1:30)
ctrl <- AddMetaData(object = ctrl, metadata = predictions)
```

```{r}
DimPlot(ctrl, group.by = "predicted.id", label = T, repel = T) + NoAxes()
```

Now plot how many cells of each celltypes can be found in each cluster.

```{r}
ggplot(ctrl@meta.data, aes(x=CCA_snn_res.0.5, fill = predicted.id)) + geom_bar() + theme_classic()
```

# scPred
scPred will train a classifier based on all principal components. First, `getFeatureSpace` will create a scPred object stored in the `@misc` slot where it extracts the PCs that best separates the different celltypes. Then `trainModel` will do the actual training for each celltype.

```{r}
reference <- getFeatureSpace(reference, "cell_type")

reference <- trainModel(reference)
```

We can then print how well the training worked for the different celltypes by printing the number of PCs used for each, the ROC value and Sensitivity/Specificity. Which celltypes do you think are harder to classify based on this dataset?

```{r}
get_scpred(reference)
```

You can optimize parameters for each dataset by chaning parameters and testing different types of models, see more at: https://powellgenomicslab.github.io/scPred/articles/introduction.html. But for now, we will continue with this model.

 Now, lets predict celltypes on our data, where scPred will align the two datasets with Harmony and then perform classification.

```{r}
ctrl <- scPredict(ctrl, reference)
```

```{r}
DimPlot(ctrl, group.by = "scpred_prediction", label = T, repel = T) + NoAxes()
```

Now plot how many	cells of each celltypes	can be found in	each cluster.

```{r}
ggplot(ctrl@meta.data, aes(x=CCA_snn_res.0.5, fill = scpred_prediction)) + geom_bar() + theme_classic()
```

# Compare results

Now we will compare the output of the two methods using the convenient function in scPred `crossTab` that prints the overlap between two metadata slots.

```{r}
crossTab(ctrl, "predicted.id", "scpred_prediction")
```


# GSEA with celltype markers

Another option, where celltype can be classified on cluster level is to use gene set enrichment among the DEGs with known markers for different celltypes. Similar to how we did functional enrichment for the DEGs in the Differential expression exercise. 
There are some resources for celltype gene sets that can be used. Such as [CellMarker](http://bio-bigdata.hrbmu.edu.cn/CellMarker/), [PanglaoDB](https://panglaodb.se/) or celltype gene sets at [MSigDB](https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).
We can also look at overlap between DEGs in a reference dataset and the dataset you are analysing. 

## DEG overlap
First, lets extract top DEGs for our Covid-19 dataset and the reference dataset.
When we run differential expression for our dataset, we want to report as many genes as possible, hence we set the cutoffs quite lenient.

```{r}
# run differential expression in our dataset, using clustering at resolution 0.5
alldata <- SetIdent(alldata,value = "CCA_snn_res.0.5")
DGE_table <- FindAllMarkers(alldata,
                               logfc.threshold = 0,
                               test.use = "wilcox",
                               min.pct = 0.1,
                               min.diff.pct = 0,
                               only.pos = TRUE,
                               max.cells.per.ident = 20,
                               return.thresh = 1,
                               assay = "RNA")

# split into a list
DGE_list <- split(DGE_table, DGE_table$cluster)

unlist(lapply(DGE_list, nrow))
```

```{r}
# Compute differential gene expression in reference dataset (that has cell annotation)
reference <- SetIdent( reference, value = "cell_type")
reference_markers <- FindAllMarkers( reference , min.pct = .1 , 
                                     min.diff.pct = .2, only.pos = T, 
                                     max.cells.per.ident = 20 ,
                                     return.thresh = 1 )

# Identify the top cell marker genes in reference dataset
# select top 50 with hihgest foldchange among top 100 signifcant genes.
reference_markers <- reference_markers [ order(reference_markers$avg_log2FC,decreasing = T), ]
reference_markers %>% 
  group_by(cluster) %>% 
  top_n(-100, p_val) %>% 
  top_n(50, avg_log2FC) -> top50_cell_selection

# Transform the markers into a list
ref_list = split(top50_cell_selection$gene, top50_cell_selection$cluster)

unlist(lapply(ref_list, length))
```

Now we can run GSEA for the DEGs from our dataset and check for enrichment of top DEGs in the reference dataset.

```{r}
suppressPackageStartupMessages(library(fgsea))

# run fgsea for each of the clusters in the list
res <- lapply(DGE_list, function(x){
  gene_rank <- setNames(x$avg_log2FC, x$gene)
  fgseaRes <- fgsea( pathways=ref_list, stats=gene_rank,nperm=10000)
  return(fgseaRes)
})
names(res) <- names(DGE_list)

# You can filter and resort the table based on ES, NES or pvalue
res <- lapply(res, function(x) {x[ x$pval < 0.1 , ]} )
res <- lapply(res, function(x) {x[ x$size > 2 , ]} )
res <- lapply(res, function(x) {x[order(x$NES,decreasing = T), ]} )
res
```

Selecing top significant overlap per cluster, we can now rename the clusters according to the predicted labels. OBS! Be aware that if you have some clusters that have non-significant p-values for all the gene sets, the cluster label will not be very reliable. Also, the gene sets you are using may not cover all the celltypes you have in your dataset and hence predictions may just be the most similar celltype.
Also, some of the clusters have very similar p-values to multiple celltypes, for instance the ncMono and cMono celltypes are equally good for some clusters.

```{r}
new.cluster.ids <- unlist(lapply(res,function(x){as.data.frame(x)[1,1]}))

alldata$ref_gsea <- new.cluster.ids[as.character(alldata@active.ident)]

cowplot::plot_grid( ncol = 2,
DimPlot(alldata,label = T,group.by = "CCA_snn_res.0.5") + NoAxes(),
DimPlot(alldata,label = T, group.by = "ref_gsea") + NoAxes())
```

Compare to results with the other celltype prediction methods in the ctrl_13 sample.

```{r, fig.width=10}
ctrl$ref_gsea = alldata$ref_gsea[alldata$orig.ident == "ctrl_13"]

cowplot::plot_grid( ncol = 3,
DimPlot(ctrl,label = T,group.by = "ref_gsea") + NoAxes() + ggtitle("GSEA"),
DimPlot(ctrl,label = T, group.by = "predicted.id") + NoAxes() + ggtitle("LabelTransfer"),
DimPlot(ctrl,label = T, group.by = "scpred_prediction") + NoAxes() + ggtitle("scPred")
)
```

## With annotated gene sets
First download celltype gene sets from CellMarker.

```{r}
# Download gene marker list
if(!dir.exists("data/CellMarker_list/")) {
  dir.create("data/CellMarker_list")
  download.file(url="http://xteam.xbio.top/CellMarker/download/Human_cell_markers.txt",
                destfile = "./data/CellMarker_list/Human_cell_markers.txt")
  download.file(url="http://xteam.xbio.top/CellMarker/download/Mouse_cell_markers.txt",
                destfile = "./data/CellMarker_list/Mouse_cell_markers.txt")
}
```

Read in the gene lists and do some filtering.

```{r}
# Load the human marker table
markers <- read.delim("data/CellMarker_list/Human_cell_markers.txt")
markers <- markers [ markers$speciesType == "Human", ]
markers <- markers [ markers$cancerType == "Normal", ]

#Filter by tissue (to reduce computational time and have tissue-specific classification)
sort(unique(markers$tissueType))
grep("blood",unique(markers$tissueType),value = T)
markers <- markers [ markers$tissueType %in% c("Blood","Venous blood",
                                                "Serum","Plasma",
                                                "Spleen","Bone marrow","Lymph node"), ]

# remove strange characters etc.
celltype_list <- lapply( unique(markers$cellName) , function(x){
  x <- paste(markers$geneSymbol[markers$cellName == x],sep=",")
  x <- gsub("[[]|[]]| |-",",",x)
  x <- unlist(strsplit( x , split = ","))
  x <- unique(x [ ! x %in% c("","NA","family") ])
  x <- casefold(x,upper = T)
})
names(celltype_list) <- unique(markers$cellName)
# celltype_list <- lapply(celltype_list , function(x) {x[1:min(length(x),50)]} )
celltype_list <- celltype_list[ unlist(lapply(celltype_list,length)) < 100 ]
celltype_list <- celltype_list[ unlist(lapply(celltype_list,length)) > 5 ]

```

```{r}
# run fgsea for each of the clusters in the list
res <- lapply(DGE_list, function(x){
  gene_rank <- setNames(x$avg_log2FC, x$gene)
  fgseaRes <- fgsea( pathways=celltype_list, stats=gene_rank,nperm=10000)
  return(fgseaRes)
})
names(res) <- names(DGE_list)

# You can filter and resort the table based on ES, NES or pvalue
res <- lapply(res, function(x) {x[ x$pval < 0.01 , ]} )
res <- lapply(res, function(x) {x[ x$size > 5 , ]} )
res <- lapply(res, function(x) {x[order(x$NES,decreasing = T), ]} )

# show top 3 for each cluster.
lapply(res,head,3)
```


#CT_GSEA8:

```{r}
new.cluster.ids <- unlist(lapply(res,function(x){as.data.frame(x)[1,1]}))
alldata$cellmarker_gsea <- new.cluster.ids[as.character(alldata@active.ident)]

cowplot::plot_grid( ncol = 2,
DimPlot(alldata,label = T,group.by = "ref_gsea") + NoAxes(),
DimPlot(alldata,label = T, group.by = "cellmarker_gsea") + NoAxes()
)
```

Do you think that the methods overlap well? Where do you see the most inconsistencies?

In this case we do not have any ground truth, and we cannot say which method performs best. You should keep in mind, that any celltype classification method is just a prediction, and you still need to use your common sense and knowledge of the biological system to judge if the results make sense.

# Save data
Finally, lets save the data with predictions.

```{r}
saveRDS(ctrl,"data/results/ctrl13_qc_dr_int_cl_celltype.rds")
```


### Session Info
***

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