cluster_naming.Rmd

---
title: "scWBM Cluster Labeling"
author: "D. Ford Hannum"
date: "5/29/2020"
output: 
        html_document:
                toc: true
                toc_depth: 3
                number_sections: false
                theme: united
                highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F)
library(Seurat)
library(ggplot2)
library(data.table)
library(SingleR)
library(MAST)
library(scRNAseq)
library(tidyr)
```

```{r reading in the data}
wbm <- readRDS('./data/wbm_clustered_filtered_named.rds')
```

# Introduction

In this document I am collecting all the information and figures that are used for naming the clusters of the Sang scWBM experiment one.

*Cluster naming workflow*

1. Using singleR to compare the cluster centroids in our analysis to reference datasets, to find the reference cell type with the highest correlation to our cell type. The datasets being compare to are: 830 microarray samples of pure mouse immune cells, generated by the Immunologic Genome Project (ImmGen, Aran et al., 2019); and 358 bulk RNA-seq samples of sorted cell populations that can be found at GEO (Benayoun et al., 2019)

2. Using singleR to compare individual cells to the reference cell types (same as step 1)

3. Using marker gene expression to either confirm or alter cluster labels generated from steps 1 and 2.


```{r singleR loading datasets, include = F, warnings = F}
m.ref <- ImmGenData()
m.ref2 <- MouseRNAseqData()
ref <- list(m.ref,m.ref2)
ref_lab <- list(m.ref$label.main, m.ref2$label.main)
```
 \newpage
 
# SingleR Cluster Naming

```{r singleR cluster labeling}
SCwbm <- as.SingleCellExperiment(wbm)

# Predicitng cluster labels

pred_cluster <- SingleR(test = SCwbm, 
                        ref = ref,
                        labels = ref_lab, 
                        method = 'cluster',
                        clusters = SCwbm$seurat_clusters)

pred_cell <- SingleR(test = SCwbm,
                     ref = ref,
                     labels = ref_lab,
                     method = 'single')
```

```{r cluster final scores data manipulation}
clst_score_final <- pred_cluster$scores

# deleting columns with all NAs
clst_score_final <- as.data.frame(clst_score_final[, colSums(is.na(clst_score_final)) != nrow(clst_score_final)])

# adding the cluster as a factor instead of just a row name
clst_score_final$cluster<- as.factor(0:12)

# labeling the colnames by the reference they came from

colnames(clst_score_final)[1:7] <- paste0(colnames(clst_score_final)[1:7], '_ref1')

colnames(clst_score_final)[8:11] <- paste0(colnames(clst_score_final)[8:11], '_ref2')

clst_score_final <- gather(clst_score_final, cell_type, score, `B cells_ref1`:`T cells_ref2`, factor_key = T)

clst_score_final$score_round <- round(clst_score_final$score,2)
```

```{r built in SingleR plot}
library(pheatmap)
plotScoreHeatmap(pred_cluster)
```

Above is a graphic from a package associated with SingleR. It is interesting but has to much going on.

```{r cluster score graph bar chart}
#levels(clst_score_final$cell_type)

clst_bar <- clst_score_final[!is.na(clst_score_final$score),]

# splitting up the cell type and reference label
clst_bar$ref <- tstrsplit(clst_bar$cell_type, '_', keep = 2)[[1]]
clst_bar$cell_type <- tstrsplit(clst_bar$cell_type, '_', keep = 1)[[1]]

clst_bar$Reference <- ifelse(clst_bar$ref == 'ref1', 'ImmGen','GEO')

ggplot(clst_bar, aes( x = cluster, y = score, fill = ref)) + 
        geom_bar(stat = 'identity', position = position_dodge()) +
        coord_flip() + 
        scale_fill_manual(values = c('darkgreen', 'orange'),
                          name = 'Reference\nDataset') + 
        theme_bw() +
        ylim(0,1) + 
        geom_text(stat = 'identity', aes(label = cell_type),
                  position = position_dodge(width = 1),
                  hjust = -.1, size = 3) + 
        xlab('Cluster') + ylab ('SingleR Final Score') +
        scale_x_discrete()
```

The final score given to each cluster. The higher score between the two references for each cluster was the final label. 

# SingleR Cell Naming

```{r sc namming}
# head(pred_cell)
# dim(wbm@meta.data)
# dim(pred_cell)
# sum(rownames(pred_cell) == rownames(wbm@meta.data))

wbm@meta.data$cell_assignment <- pred_cell$pruned.labels

df <- as.data.frame(table(wbm@meta.data$seurat_clusters, wbm@meta.data$cell_assignment))

colnames(df) <- c('Cluster','Assigned Cell Type','Count')

cluster_counts <- c(summary(as.factor(wbm@meta.data$seurat_clusters)))
cluster_counts
dim(df)[1]/length(cluster_counts)
df$clst_counts <- rep(cluster_counts,18)
df$cell_perc <- round(df$Count/df$clst_counts,2)*100

new_levels <- levels(as.factor(df$`Assigned Cell Type`))[c(8,1,3,12,10,7,2,17,14,16,13,9,4,6,5,11,15,18)]
df$cell_type <- factor(df$`Assigned Cell Type`, levels = new_levels)

ggplot(df, aes(x = cell_type, y = Cluster, fill = cell_perc)) + 
        geom_tile() +
        scale_fill_gradient2(high = 'darkred', low = 'white', mid = 'red', midpoint = 50,
                             limit = c(0,100), space ="Lab",
                             name = 'Percentage of Cells\n') +
        ylab('Cluster Labels') + xlab('Cell Labels') +
        ggtitle ("Heatmap of Individual Cell Labels within Cluster Labels") + 
        theme (plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 45, hjust = 1)) +
        coord_fixed()
```

# Lineage Specific Markers

Next I'm going to use lineage specific markers to manually look at the cell clusters. List was provided by Dr. Sang.

* **Granulocytes**: *Elane, Mpo, Ctsg, Prtn3, Azu1*

* **B-cell**: *Ighd, Cd20, Cd22, Jchain*

* **Megakaryocytes**: *Itga2b, Gp9, Pf4, Selp, Gp1ba*

* **HSPC**: *Crhbp, Emcn, Hlf, Avp, Cd34, Kit, Sca-1*

* **Monocyte**: *Cd14*

* **Macrophage**: *Cd45, F480, Lyz2*

* **Erythroid**: *Epor, Klf1, Tfr2, Csf2rb, Gypa*

* **T-cell/NK**: *Cd3g*

* **MEP**: *Lbg, Nr4a1, Gpr141, Gata1*

## Granulocytes

```{r changing wbm identi}
# changing the identity to the cluster number again
new_cluster_ids <- 0:12
names(new_cluster_ids) <- levels(wbm)
wbm <- RenameIdents(wbm, new_cluster_ids)
```

```{r gran markers}
gran_markers <- c('Elane', 'Mpo', 'Ctsg', 'Prtn3', 'Azu1')
gran_markers %in% rownames(wbm)
gran_markers <- gran_markers[gran_markers %in% rownames(wbm)]

VlnPlot(wbm, features = gran_markers,
        pt.size = 0)

FeaturePlot(wbm, features = gran_markers)
```

These markers are only consistently expressed in cluster 6. Some marker expression in clusters 4, 5 and 7.

\newpage

## B-cell

```{r b cell markers}
b_markers <- c('Ighd', 'Cd20', 'Cd22', 'Jchain')
b_markers <- b_markers[b_markers %in% rownames(wbm)]

VlnPlot(wbm, features = b_markers,
        pt.size = 0)

FeaturePlot(wbm, features = b_markers)
```
Jchain showed little to no expression in any cluster.

Expression in clusters 3 and 10, almost exclusively.

\newpage 

## Megakaryocytes

```{r mk markers}
mk_markers <- c('Itga2b', 'Gp9', 'Pf4', 'Selp', 'Gp1ba')
mk_markers <- mk_markers[mk_markers %in% rownames(wbm)]

VlnPlot(wbm, features = mk_markers,
        pt.size = 0)

FeaturePlot(wbm, features = mk_markers)
```

Most expression was found in cluster 12, with strong Itga2b expression in 4 and some Pf4 expression in 8.

\newpage

## HSPC

```{r hspc markers}
hspc_markers <- c('Crhbp', 'Emcn', 'Hlf', 'Avp', 'Cd34', 'Kit', 'Sca-1')
hspc_markers <- hspc_markers[hspc_markers %in% rownames(wbm)]

VlnPlot(wbm, features = hspc_markers, pt.size = 0)
```
Emcn, Avp, Cd34, and Hlf showed little to no expression in any cluster.

Crhbp has mild expression in cluster 10, and Kit has mild expression in cluster 6.

\newpage

## Monocyte

```{r mono marker}
mono_marker <- 'Cd14'
#'Cd14' %in% rownames(wbm)

VlnPlot(wbm, features = mono_marker, pt.size = 0) + NoLegend()
```

Cd14 showed little to no expression

\newpage 

## Macrophage

```{r macro markers}
macro_markers <- c('Cd45', 'F480', 'Lyz2')
macro_markers <- macro_markers[macro_markers %in% rownames(wbm)]

VlnPlot(wbm, features = macro_markers, pt.size = 0)
```

Lyz2 is widely expressed in all clusters

\newpage

## Erythroid

```{r ery markers}
ery_markers <- c('Epor', 'Klf1', 'Tfr2', 'Csf2rb', 'Gypa')
ery_markers <- ery_markers[ery_markers %in% rownames(wbm)]

VlnPlot(wbm, features = ery_markers, pt.size = 0)
FeaturePlot(wbm, features = ery_markers)
```
Epor, Klf1, and Tfr2 showed little to no expression in any clusters.

Csf2rb shows widespread expression, with highest expression in cluster 4. Gypa only shows expression in cluster 9

\newpage

## T-cell/NK

```{r tnk markers}
tnk_markers <- c('Cd3g')
tnk_markers <- tnk_markers[tnk_markers %in% rownames(wbm)]

VlnPlot(wbm, features = tnk_markers, pt.size =0)
FeaturePlot(wbm, features = tnk_markers)
```

Cd3g shows expression in cluster 11, almost exclusively.

\newpage

## MEP

```{r mep markers}
mep_markers <- c('Lbg', 'Nr4a1', 'Gpr141', 'Gata1')
mep_markers <- mep_markers[mep_markers %in% rownames(wbm)]

VlnPlot(wbm, features = mep_markers, pt.size = 0)
```

Gpr141 shows expression in multiple clusters. Expression of Gata1 is exclusive to cluster 4.

# Lineage Specfic Markers cont.

These are markers that I found through literature search and though looking at [pangloaddb](pangloadb.se) (for example top markers for [HSPCs](https://panglaodb.se/markers.html?cell_type=%27Hematopoietic%20stem%20cells%27))

Many of the markers are the same between my list and Dr. Sang's list.

```{r reading in marker list}
markers <- read.csv('./marker_list.csv', header = T)
markers
```

\newpage

## B-cells
```{r otro b-cell markers}
mrkrs <- markers$B.cell
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

Clusters 3 and 10 are high in B-cell markers which is consistent through all the stages

\newpage

## MKs
```{r otro mk markers}
mrkrs <- markers$MK
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## HSPCs
```{r otro hspc markers}
mrkrs <- markers$HSPC
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```


\newpage

## Monocytes
```{r otro mono markers}
mrkrs <- markers$Monocyte
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```


\newpage

## Macrophages
```{r otro macrophage markers}
mrkrs <- markers$Macrophage
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## Erythroid
```{r otro ery markers}
mrkrs <- markers$Erythroid
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## T-cell/NK
```{r otro tnk markers}
mrkrs <- markers$T.cell.NK
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## MEP
```{r otro mep markers}
mrkrs <- markers$MEP
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## Proliferation 
```{r otro pro markers}
mrkrs <- markers$ProliferationMarkers
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## Myeloid
```{r otro myeloid markers}
mrkrs <- markers$Myeloid
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```


\newpage

## Erythroid Prog.
```{r otro ery pro markers}
mrkrs <- markers$Eprog
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}
```

\newpage

## Plasma
```{r otro plasma markers}
mrkrs <- markers$Plasma
mrkrs <- mrkrs[mrkrs %in% rownames(wbm)]

VlnPlot(wbm, features = mrkrs, pt.size = 0)

loops <- ceiling(length(mrkrs) / 4)
for (i in 1:loops){
        start <- 4*(i-1) + 1
        end <- 4*i
        
        if (end > length(mrkrs)){
                mrkrs2 <- mrkrs[start:length(mrkrs)]
        }
        else{
                mrkrs2 <- mrkrs[start:end]
        }
        
        print(FeaturePlot(wbm, features = mrkrs2, ncol = 2))
}

```