---
title: "Supereffectors_PART3"
author: "Veronika Niederlova"
date: "15 10 2021"
output: 
  html_document: 
    fig_height: 2.8
    fig_width: 4
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(include = TRUE, warning = FALSE, message = FALSE, error = TRUE, cache = TRUE)

library("Matrix")
library(Seurat)
library(DT)
library(dplyr)
library(here)
library(ggplot2)
library(kableExtra)
library(cowplot)
library(tidyverse)
library(reshape)
library(SingleR)
library(ProjecTILs)
library(ggprism)
library(cowplot)
library(fgsea)
library(org.Mm.eg.db)

```

The input of this analysis is a processed Seurat file output of Supereffectors_Analysis_PART2. The file is called exp07_W3_Treg_processed_data_set_2_filtered.rds and we expect that it's located in the ./data folder. In this part of the analysis workflow, we will share the code that was used for generating all figures in the manuscript.

## Load the data
```{r}
path_to_downloaded_file <- "G://48_lab/Project Tregs DIA/! Manuscript in preparation/RNAseq DATA for upload/exp07_W3_Treg_processed_data_set_2_filtered.rds"

seu_wocc <- readRDS(path_to_downloaded_file)
```


## Figure 5

### Figure 5B

•	UMAP plot with updated clustering 

```{r fig.width=3.3}

DimPlot(seu_wocc, group.by = "seurat_clusters2", pt.size = 0.5) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + ggtitle("")

# ggsave(filename = "final_fig/fig_5b_leg.png", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5b_leg.eps", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5b_leg.svg", width = 16, height = 13, units = "cm")



```

### Figure 5C

•	UMAP with KLRK1 expression 

```{r fig.width=3.3}
FeaturePlot(seu_wocc, features = "Klrk1", pt.size = 0.5) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + scale_x_continuous(breaks = c(-6,-3,0,3,6))

# ggsave(filename = "final_fig/fig_5c_leg.png", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5c_leg.eps", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5c_leg.svg", width = 16, height = 13, units = "cm")

```

### Figure 5D

•	UMAP with differentiated cells from DEREG- and DEREG+ 

```{r, fig.width=3.7}
Idents(seu_wocc) <- "dereg"
DimPlot(seu_wocc,  pt.size = 0.5, cols = c("indianred2","royalblue")) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + ggtitle("")
# ggsave(filename = "final_fig/fig_5d_leg.png", width = 18.3, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5d_leg.eps", width = 18.3, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5d_leg.svg", width = 18.3, height = 13, units = "cm")
```



### Figure 5E

•	Plots showing frequency of cells of different clusters in DEREG+ and DEREG- samples 

```{r fig.width=9}

# Calculate frequencies of clusters in different anuimals
df4 <- seu_wocc@meta.data %>% group_by(hashtags, seurat_clusters2) %>% 
  summarise(n = n()) %>%
  mutate(freq = n / sum(n)) %>% 
  mutate(dereg = case_when(hashtags %in% c("Mouse_Hashtag_5","Mouse_Hashtag_6","Mouse_Hashtag_7") ~ "DEREG-",
                                                                  TRUE ~ "DEREG+")) 
# Add p-values
labels <- c()
for(i in 0:4){
  df <- df4 %>% filter(seurat_clusters2 == i)
  pval <- t.test(df$freq~df$dereg)$p.value*5
  
  if(pval>1){
    pval <- "1"
  } else {if(pval<0.001){
    pval <- "<0.001"
  } else { pval <-round(pval, digits = 3)}}
  
  label_new <- paste("Cluster",i,"\n",pval)
  labels <- c(labels,label_new)
}

# Plot results
 df4 %>%
   ggplot(aes(x = dereg, y = freq*100)) + 
 geom_dotplot(binaxis='y', stackdir='center', dotsize=0) + 
   geom_jitter(position=position_jitter(0.2), size = 3, aes(color = dereg)) +
   theme_minimal() + 
   facet_wrap(~factor(seurat_clusters2, labels = labels), scales = "free", ncol = 5) + 
   ylim(0,NA) + 
   theme_prism(base_fontface = "plain", base_line_size = 0.8) + 
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11), 
         axis.title = element_text(size = 11), 
         legend.text = element_text(size = 11)) + 
   xlab("") + scale_color_manual(values = c("royalblue3","indianred2")) +
   stat_summary(fun.y = median, fun.ymin = median, fun.ymax = median,
                 geom = "crossbar", width = 1) + theme(axis.text.x = element_blank())
 
# Save figures
# ggsave(filename = "final_fig/fig_5e_leg.png", width = 38, height = 10, units = "cm")
# ggsave(filename = "final_fig/fig_5e_leg.svg", width = 38, height = 10, units = "cm")

```

### Figure 5F

•	reclustering of cluster 4 (SE cluster)

```{r}
se <- subset(seu_wocc, seurat_clusters2 == 4)

se <- ScaleData(se, verbose = FALSE, vars.to.regress = c("S.Score", "G2M.Score"))
se <- FindVariableFeatures(se, nfeatures = 100, verbose = FALSE)
se <- RunPCA(se, npcs = 7, verbose = FALSE)
se <- RunUMAP(se, reduction = "pca", dims = 1:7)
se <- FindNeighbors(se, dims = 1:7)
se <- FindClusters(se, resolution = 0.2, n.start = 20, n.iter = 50)

DimPlot(se, cols = c("orange2","forestgreen")) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + 
  ggtitle("") + scale_x_continuous(breaks = c(-5,-2.5,0,2.5,5))+xlim(c(-6,6))


# ggsave(filename = "final_fig/fig_5f1_leg.png", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f1_leg.eps", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f1_leg.svg", width = 16, height = 13, units = "cm")

DimPlot(se, cols = c("royalblue3","indianred2"), group.by = "dereg") + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + 
  ggtitle("") + scale_x_continuous(breaks = c(-5,-2.5,0,2.5,5)) +xlim(c(-6,6))
# ggsave(filename = "final_fig/fig_5f2_leg.png", width = 18.3, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f2_leg.eps", width = 18.3, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f2_leg.svg", width = 18.3, height = 13, units = "cm")
```

### Figure 5G

•	SE subcluster, expression of Itga4 and Il7r

```{r}
FeaturePlot(se, features = "Itga4", max.cutoff = 2) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11))  + xlim(c(-6,6))
# ggsave(filename = "final_fig/fig_5f3_leg.png", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f3_leg.eps", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f3_leg.svg", width = 16, height = 13, units = "cm")

FeaturePlot(se, features = "Il7r", max.cutoff = 2) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + xlim(c(-6,6))
# ggsave(filename = "final_fig/fig_5f4_leg.png", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f4_leg.eps", width = 16, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_5f4_leg.svg", width = 16, height = 13, units = "cm")


```



## Figure 5H

• Projection of the supereffector subclusters onto the original UMAP plot.


```{r}
md <- seu_wocc@meta.data
md$barcode <- colnames(seu_wocc)
md <- md %>% mutate(is_subcluster = case_when(barcode %in% colnames(subset(se,seurat_clusters ==1)) ~ "Subcluster 1",
                                              barcode %in% colnames(subset(se,seurat_clusters ==0)) ~ "Subcluster 0",
                                              TRUE ~ "Other"))

seu_wocc@meta.data <- md

DimPlot(seu_wocc, cells.highlight = list(colnames(subset(se,seurat_clusters ==0)),colnames(subset(se,seurat_clusters ==1))), cols.highlight = c("forestgreen","orange2")) + theme(axis.text = element_text(size = 11), axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + ggtitle("")
# ggsave(filename = "final_fig/fig_s5e5_leg.png", width = 18.3, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_s5e5_leg.eps", width = 18.3, height = 13, units = "cm")
# ggsave(filename = "final_fig/fig_s5e5_leg.svg", width = 18.3, height = 13, units = "cm")

```


# Figure S5


### Figure S5B

•	Expression of GzmB in DEREG- and DEREG+ samples (Fig. S5 B), change colours – red for DEREG+, blue for DEREG- 

```{r}
VlnPlot(seu_wocc, features = "Gzmb", cols = c("indianred2","royalblue3"), flip = TRUE) + theme(axis.text.x = element_text(angle = 0, size = 11), axis.text = element_text(size = 11), axis.title = element_text(size = 11), legend.text = element_text(size = 11))

# ggsave(filename = "final_fig/fig_s5b_leg.png", width = 16, height = 8, units = "cm")
# ggsave(filename = "final_fig/fig_s5b_leg.eps", width = 16, height = 8, units = "cm")
# ggsave(filename = "final_fig/fig_s5b_leg.svg", width = 16, height = 8, units = "cm")
```

### Figure S5C

•	Expression of different markers 

```{r fig.height=13, fig.width=28}
FeaturePlot(seu_wocc, features = c("Gzma", "Gzmb", "Gzmk", "Il2ra", "Il7r", "Cd7", "Ifitm1", "Ifitm2", "Ifitm3", "Klrd1", "Id2", "Id3", "Itga4", "Itgae", "Sell", "Ccr7", "Cxcr6", "Cd27"), ncol = 6) + theme(axis.text.x = element_text(angle = 0, size = 16), axis.text = element_text(size = 11), axis.title = element_text(size = 11), legend.text = element_text(size = 11))

# ggsave(filename = "final_fig/fig_s5b_leg.png", width = 16, height = 8, units = "cm")
# ggsave(filename = "final_fig/fig_s5b_leg.eps", width = 16, height = 8, units = "cm")
# ggsave(filename = "final_fig/fig_s5b_leg.svg", width = 16, height = 8, units = "cm")
```


### Figure S5D

•	ProjecTILs projection of supereffector cluster onto LCMV data. See the original publication [Andreatta et al., Nat Comm, 2021](https://www.nature.com/articles/s41467-021-23324-4). Classifications are stored in the object, but we provide full code for their calculation.

```{r fig.height=16, fig.width=28, cache.lazy=FALSE}

Tregs_filt_wocc <- readRDS("G://48_lab/Project Tregs DIA/! Manuscript in preparation/RNAseq DATA for upload/exp07_W3_Treg_processed_data_set_2_filtered.rds")

Idents(Tregs_filt_wocc) <- Tregs_filt_wocc$seurat_clusters2

ref <- readRDS("G://48_lab/Project scRNAseq/Scripts Verca/Supereffectors_scRNAseq/data/ref_LCMV_Atlas_mouse_v1.rds")

query0 <- subset(Tregs_filt_wocc, seurat_clusters2 == 0)
query.projected0 <- make.projection(query0, ref=ref, skip.normalize = T)
p0 <- plot.projection(ref, query.projected0, linesize = 0.6) +
  ggtitle("Cluster 0") +
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 11),
         legend.text = element_text(size = 11))

query1 <- subset(Tregs_filt_wocc, seurat_clusters2 == 1)
query.projected1 <- make.projection(query1, ref=ref, skip.normalize = T)
p1 <- plot.projection(ref, query.projected1, linesize = 0.6) +
  ggtitle("Cluster 1") +
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 11),
         legend.text = element_text(size = 11))

query2 <- subset(Tregs_filt_wocc, seurat_clusters2 == 2)
query.projected2 <- make.projection(query2, ref=ref, skip.normalize = T)
p2 <- plot.projection(ref, query.projected2, linesize = 0.6) +
  ggtitle("Cluster 2") +
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 11),
         legend.text = element_text(size = 11))

query3 <- subset(Tregs_filt_wocc, seurat_clusters2 == 3)
query.projected3 <- make.projection(query3, ref=ref, skip.normalize = T)
p3 <- plot.projection(ref, query.projected3, linesize = 0.6) +
  ggtitle("Cluster 3") +
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 11),
         legend.text = element_text(size = 11))

query4 <- subset(Tregs_filt_wocc, seurat_clusters2 == 4)
 query.projected4 <- make.projection(query4, ref=ref, skip.normalize = T)
p4 <- plot.projection(ref, query.projected4, linesize = 0.6) +
  ggtitle("Cluster 4") +
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 11),
         legend.text = element_text(size = 11))

cowplot::plot_grid(p0,p1,p2,p3,p4, ncol = 5)

# ggsave(filename = "final_fig/fig_s5d_leg.png", width = 74, height = 20, units = "cm")
# ggsave(filename = "final_fig/fig_s5d_leg.svg", width = 74, height = 20, units = "cm")
```


## Figure S5F

•	Expression of Klrk1 and Cd49d in one plot (flow gating)

```{r fig.height=4, fig.width=14}

FeaturePlot(seu_wocc,  pt.size = 0.5, features = c("Klrk1","Itga4"), blend = T, cols = c("indianred2","blue2"), max.cutoff = c(1,2), blend.threshold = 0) + theme(axis.text = element_text(size = 11), 
                                                                        axis.title = element_text(size = 11),
                                                                        legend.text = element_text(size = 11)) + ggtitle("")

```

Fig S5 - not used

•	SE subcluster, percentage in DEREG+ and DEREG- animals

```{r fig.width=6.5}
df5 <- se@meta.data %>% group_by(hashtags, seurat_clusters) %>%
  summarise(n = n()) %>%
  mutate(freq = n / sum(n)) %>%
  mutate(dereg = case_when(hashtags %in% c("Mouse_Hashtag_5","Mouse_Hashtag_6","Mouse_Hashtag_7") ~ "DEREG-",
                                                                  TRUE ~ "DEREG+"))

 df5 %>%
   ggplot(aes(x = dereg, y = freq*100)) +
 geom_dotplot(binaxis='y', stackdir='center', dotsize=0) +
   geom_jitter(position=position_jitter(0.2), size = 3, aes(color = dereg)) +
   theme_minimal() +
   facet_wrap(~factor(seurat_clusters, labels = c("Subcluster 0", "Subcluster 1")), scales = "free", ncol = 5) +
   ylim(0,NA) +
   theme_prism(base_fontface = "plain", base_line_size = 0.8) +
   theme(axis.text.x = element_text(size = 11),
         title = element_text(size = 11),
         axis.text = element_text(size = 11),
         axis.title = element_text(size = 11),
         legend.text = element_text(size = 11)) +
   xlab("") + scale_color_manual(values = c("royalblue3","indianred2")) +
   stat_summary(fun.y = median,
                 geom = "crossbar", width = 1)
 
 
# ggsave(filename = "final_fig/fig_5g_leg.png", width = 16, height = 7, units = "cm")
# ggsave(filename = "final_fig/fig_5g_leg.svg", width = 16, height = 7, units = "cm")
```

### SessionInfo

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