---
title: "Marker Identification & Cluster Annotation Helpers"
date: 'Compiled: `r format(Sys.Date(), "%B %d, %Y")`'
output: rmarkdown::html_vignette
theme: united
df_print: kable
vignette: >
  %\VignetteIndexEntry{Marker Identification & Cluster Annotation Helpers}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
***

<style>
p.caption {
  font-size: 0.9em;
}
</style>
 
```{r setup, include=FALSE}
all_times <- list()  # store the time for each chunk
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      now <<- Sys.time()
    } else {
      res <- difftime(Sys.time(), now, units = "secs")
      all_times[[options$label]] <<- res
    }
  }
}))
knitr::opts_chunk$set(
  tidy = TRUE,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  time_it = TRUE
)
```



# Marker Identification & Cluster Annotation  
scCustomize has several helper functions to assist with identification of marker genes and annotation of clusters.
Load packages & Data
```{r init}
# Load Packages
library(ggplot2)
library(dplyr)
library(magrittr)
library(Seurat)
library(scCustomize)
library(qs)

# Load example dataset for tutorial
pbmc <- pbmc3k.SeuratData::pbmc3k.final
```

## Helper functions available are:

* `Extract_Top_Markers()`: Extract either a filtered data.frame or vector of gene IDs (with option to export named vector) containing top N markers from data.frame.
* `Create_Cluster_Annotation_File()`: Create simple csv file for use in annotating clusters.
* `Pull_Cluster_Annotation()`: Pull cluster annotation information quickly from data.frame/tibble or annotation csv file.
* `Rename_Clusters()`: Wrapper function to easily rename clusters.


We will use pbmc3K dataset from SeuratData package as example for this vignette.  

```{r include=FALSE}
# Update pbmc check
pbmc <- UpdateSeuratObject(pbmc)
```

## Percent Difference in Expression
If we run `FindAllMarkers()` like normal, this is the output.
```{r}
# Run basic FindAllMarkers DE test
all_markers <- FindAllMarkers(object = pbmc, verbose = F)

all_markers <- all_markers %>%
  group_by(cluster) %>%  
  arrange(desc(avg_log2FC)) %>%   
  arrange(cluster)
```


```{r echo=FALSE}
head(all_markers, 5) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 
```

While `Seurat::FindAllMarkers()`returns the percent of cells in identity 1 (pct.1) and identity 2 (pct.2) that express a marker it can be helpful to view the difference in these two measures in addition to the values alone.  
  
scCustomize contains helper function: `Add_Pct_Diff()` to add the percent difference between two clusters.  `Add_Pct_Diff` can be used with any output from Seurat, SeuratWrappers, or other DE tests and can be piped with the function itself.

### Use `Add_Pct_Diff` piped with the marker function.
```{r}
# run piped with Seurat command
all_markers_pct <- FindAllMarkers(object = pbmc, verbose = F) %>% 
  Add_Pct_Diff()

all_markers_pct <- all_markers_pct %>% 
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>% 
  arrange(cluster)
```

```{r echo=FALSE}
head(all_markers_pct, 5) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 

```


### Use `Add_Pct_Diff` post-hoc
`Add_Pct_Diff` can also simply be run after generating marker data.frame:
```{r eval=FALSE}
all_markers_pct <- Add_Pct_Diff(marker_dataframe = all_markers)
```

### Use with non-Seurat data.frames
`Add_Pct_Diff` can be used with any data.frame as long as the column names that specify percent expressed in the 2 groups are provided
```{r eval=FALSE}
all_markers_pct <- Add_Pct_Diff(marker_dataframe = all_markers, pct.1_name = "percent1", pct.2_name = "percent2")
```



## Extract the top N marker genes per cluster, identity, treatment etc
scCustomize contains function `Extract_Top_Markers` to get the top N genes from each cluster (or other group of choice). It has the option to return either a vector of genes or a filtered data.frame.

### Returns a Named Vector (Default)
```{r}
# By default it returns a named vector of genes (each gene named for cluster it was associated with)
top_5 <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, rank_by = "avg_log2FC")
head(top_5, 10) 
```

### Return data.frame (Optional)
```{r eval=FALSE}
top_5_df <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, data_frame = TRUE, rank_by = "avg_log2FC")

head(top_5_df, 10) 
```

```{r echo=FALSE}
top_5_df <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, data_frame = TRUE, rank_by = "avg_log2FC")

head(top_5_df, 10) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 
```


### Use with non-Seurat data.frames
`Extract_Top_Markers` uses the data.frame output of `Seurat::FindAllMarkers()` by default but can be used with any data.frame by supplying few extra parameters.  Let's use the function on the output from [presto package](https://github.com/immunogenomics/presto) and order the results by auc instead of log2FC.
```{r}
# Read in pre-computed data.frame results
presto_markers <- qread("assets/presto_markers.qs")

# Extract top markers
top_5_presto <- Extract_Top_Markers(marker_dataframe = presto_markers, num_genes = 5, group_by = "cluster", rank_by = "auc", gene_column = "gene")
head(top_5_presto, 10)
```

### Unique Vectors
When returning an unnamed vector you may want to return a unique list without any duplicate features.  This is often true when supplying a list of features to plotting function like `DoHeatmap()` or `Cluster_DotPlot()`.  You can do this by setting parameter `make_unique = TRUE`.

```{r}
# Default parameters
top_5 <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, rank_by = "avg_log2FC", named_vector = FALSE)

any(duplicated(x = top_5))

# Set `make_unique = TRUE`
top_5_unique <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, rank_by = "avg_log2FC", named_vector = FALSE, make_unique = TRUE)

any(duplicated(x = top_5_unique))
```


## Pairing `Extract_Top_Markers` with `Iterate_FeaturePlot_scCustom()`
The advantage of exporting a named vector from `Extract_Top_Markers` is that `Iterate_FeaturePlot_scCustom` can detect the names and either add those to the title of the plot (if `single_pdf = TRUE`) or the file names if `FALSE`.
```{r eval=FALSE}
Iterate_FeaturePlot_scCustom(seurat_object = pbmc, gene_list = top_5_presto, single_pdf = T)
```

<center>
<img src="../reference/figures/assets/Iterate_named_plots.png" width=49%>  

*Example plot showing the gene name and additional name of the plot.*
</center>

## Create simple annotation files
Now for a simple helper we can use `Create_Cluster_Annotation_File` to create simple csv file to use for annotation.
```{r eval=FALSE}
Create_Cluster_Annotation_File(file_path = "data/", file_name = "ExpXX_Round01_Annotation")
```

Here is an example that has already been filled out (using Excel, Numbers, or other editor program).
```{r echo=FALSE}
annotation_example <- read.csv(file = "assets/ExpXX_Round01_Annotation_Completed.csv", stringsAsFactors = FALSE)
```

```{r echo=FALSE}
head(annotation_example, 5) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 
```


## Extract annotation information
There is also complementary function to pull information back out of the this (or other) csv files or data.frames regarding clustering annotation.
```{r}
annotation_info <- Pull_Cluster_Annotation(annotation = "assets/ExpXX_Round01_Annotation_Completed.csv")
```
By default the function expects data.frame or file in same format as `Create_Cluster_Annotation_File` but has optional parameters that can be adjusted to file/data.frame with different column headers (see `cluster_name_col` and `cell_type_col`).

### `Pull_Cluster_Annotation` results
The result of `Pull_Cluster_Annotation` is a list that contains:  

* One entry for each unique cell type.  The entry is named for the cell type and is a vector of the numbered clustered that correspond to that cell type.
* One entry called "new_cluster_idents" containing a vector of idents in order that can be used for renaming clusters either using `Rename_Clusters` function (see below) or manually.
* One entry containing vector of cell sub-types (if present).

<center>
<img src="../reference/figures/assets/annotation_info.png">
*Example output of `Pull_Cluster_Annotation*
</center>

## Renaming Clusters
Finally can use `Rename_Clusters` to easily rename and set active.ident in Seurat Object.  This is simple wrapper around Seurat commands but adds additional checks/warnings and is simple one-liner.
```{r eval=FALSE}
obj_renamed <- Rename_Clusters(seurat_object = obj, new_idents = annotation_info$new_cluster_idents)
```
`Rename_Clusters` can take input from either `Pull_Cluster_Annotation` "new_cluster_idents" or any correctly ordered vector of new idents.  

* If vector is not yet named (with the current levels of the Seurat Object) then `Rename_Clusters` will perform that step.  If already named this step will be omitted.
* `Rename_Clusters` contains optional parameter to create new column in meta.data slot to store the old idents.  Simply supply desired meta.data column name to `meta_col_name` parameter.


