---
title: "Sintonetal_Full Script"
output: html_document
date: '2023-05-15'
---
This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. 

### COMPLETE ANALYSIS OF MICE HYPOTHALAMUS SCRNASEQ ###

#install.packages("Seurat")
#BiocManager::install("scater")
#BiocManager::install("biomaRt")
#install.packages("tidyverse")
#install.packages('metap')
#install.packages("dplyr")

```{r}
##Load packages
library(Seurat)
library(scater)
library(biomaRt)
library(tidyverse)
library(dplyr)
library(patchwork)
library(metap)
```

##Seurat
```{r}
##Load Seurat objects
load("Naive")
load("Infected")

##Add groups
Naive$Group <- "Naive"
Infected$Group <- "Infected"

##Merge dataset
merged <- merge(Naive, c(Infected), 
                add.cell.ids = c("Naive", "Infected"))
```

```{r}
##Calculate UMI, mitochondrial/ribosomal and T. brucei transcripts percentage
merged$log10GenesPerUMI <- log(merged$nFeature_RNA)/log(merged$nCount_RNA)
merged <- PercentageFeatureSet(merged, pattern = "^mt-", col.name = 'percent.mt')
merged <- PercentageFeatureSet(merged, pattern = "^Rp[sl]", col.name = 'percent.ribo')
merged <- PercentageFeatureSet(merged, pattern = "^Tb[0-9]", col.name = "t.brucei") ##No parasite genes found

##Visualise as violin plots
plot1.merged <- FeatureScatter(object = merged, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "orig.ident")
plot2.merged <- FeatureScatter(object = merged, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "orig.ident")
plot3.merged <- FeatureScatter(object = merged, feature1 = "nCount_RNA", feature2 = "percent.ribo", group.by = "orig.ident")
plot1.merged + plot2.merged + plot3.merged + plot_layout(guides = 'collect')

plot4.merged <- VlnPlot(merged, features = "nFeature_RNA", pt.size = 0.5, group.by = "orig.ident") +
  scale_fill_manual(values = wes_palette(6, name = "FantasticFox1", type = "continuous"), name = "")
plot5.merged <- VlnPlot(merged, features = "nCount_RNA", pt.size = 0.5, group.by = "orig.ident") +
  scale_fill_manual(values = wes_palette(6, name = "FantasticFox1", type = "continuous"), name = "")
plot6.merged <- VlnPlot(merged, features = "percent.mt", pt.size = 0.5, group.by = "orig.ident") +
  scale_fill_manual(values = wes_palette(6, name = "FantasticFox1", type = "continuous"), name = "")
plot7.merged <- VlnPlot(merged, features = "percent.ribo", pt.size = 0.5, group.by = "orig.ident") +
  scale_fill_manual(values = wes_palette(6, name = "FantasticFox1", type = "continuous"), name = "") 
plot8.merged <- VlnPlot(merged, features = "log10GenesPerUMI", pt.size = 0.5, group.by = "orig.ident") +
  scale_fill_manual(values = wes_palette(6, name = "FantasticFox1", type = "continuous"), name = "") 

plot4.merged + NoLegend() + labs(x = "Group", y = "Number of Genes") + theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x.bottom = element_text(size =20), axis.title.y.left = element_text(size =20))
plot5.merged + NoLegend() + labs(x = "Group", y = "Number of UMIs") + theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x.bottom = element_text(size =20), axis.title.y.left = element_text(size =20))
plot6.merged + NoLegend() + labs(x = "Group", y = "% of mitochondrial gene transcripts") + theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x.bottom = element_text(size =20), axis.title.y.left = element_text(size =20))
plot7.merged + NoLegend() + labs(x = "Group", y = "% of ribosomal gene transcripts") + theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x.bottom = element_text(size =20), axis.title.y.left = element_text(size =20))
plot8.merged + NoLegend() + labs(x = "Group", y = "Log10 Genes Per UMI") + theme(axis.text.x = element_text(size = 20), axis.text.y = element_text(size = 20), axis.title.x.bottom = element_text(size =20), axis.title.y.left = element_text(size =20))
```

```{r}
##Filter the data to remove cells with unusual QC metrics
Filtered <- subset(merged, subset = nFeature_RNA > 200 & nFeature_RNA < 2000
                   & nCount_RNA > 100 & nCount_RNA < 20000 & percent.mt < 5 & 
                     percent.ribo < 45)
```

```{r}
##Normalisation and scaling
Filtered <- NormalizeData(Filtered)
Filtered <- FindVariableFeatures(Filtered, selection.method = "vst")
Filtered <- ScaleData(Filtered, features = rownames(Filtered))
Filtered <- RunPCA(Filtered, features = VariableFeatures(Filtered), nfeatures.print = 10)
```

```{r}
##Cell cycle scoring
##Convert Human Gene list to mus musculus
convertHumanGeneListtoMM <- function(x){require("biomaRt")
  human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
  mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
  
  genesV2 = getLDS(attributes = c("hgnc_symbol"), filters = "hgnc_symbol", values = x , mart = human, attributesL = c("mgi_symbol"), martL = mouse, uniqueRows=T)
  
  humanx <- unique(genesV2[, 2])
  
  # Print the first 6 genes found to the screen
  print(head(humanx))
  return(humanx)
}

##Extract S phase markers
m.s.genes <- convertHumanGeneListtoMM(cc.genes.updated.2019$s.genes)
##Extract G2/M phase markers
m.g2m.genes <- convertHumanGeneListtoMM(cc.genes.updated.2019$g2m.genes)
##Perform cell cycle scoring
Filtered <- CellCycleScoring(object = Filtered, g2m.features = m.g2m.genes, s.features = m.s.genes, set.ident = TRUE)
```

##Scater
```{r}
##Perform QC for cells and genes
filt.sce <- as.SingleCellExperiment(Filtered) ##Convert Seurat object to sce 
filt.sce <- addPerCellQC(filt.sce, subsets=list(Mito=grep("^mt-", rownames(filt.sce)))) ##Cell QC
filt.sce <- addPerFeatureQC(filt.sce) ##Genes QC
plotColData(filt.sce, x = "sum", y="detected", colour_by= "Sample") ##Visualise as a scatter plot

##Top 25 most expressed features
filt.sce.HE <- plotHighestExprs(filt.sce, exprs_values = "counts", n = 25) +
  theme(text = element_text(size=6))
```

##Seurat - continuation
```{r}
##Removal of overexpressed gene and parasite genes
Filtered <- Filtered[!grepl("Malat1", rownames(Filtered)), ]
Filtered <- Filtered[!grepl("Gm42418", rownames(Filtered)), ]
```

```{r}
##Split the merged Seurat object to run normalisation for each sample
Filtered.list <- SplitObject(Filtered, split.by = "orig.ident")

##SCTransform - normalisation
sctransform <- function(X){
  SCTransform(X, vars.to.regress = c('percent.mt', 'percent.ribo', 'S.Score', 'G2M.Score', 'nCount_RNA'), 
              method = "glmGamPoi")}
Filtered.list <- lapply(X = Filtered.list, FUN = sctransform) 
pca.for.dim <- function(X){RunPCA(X)}
Filtered.list <- lapply(X = Filtered.list, FUN = pca.for.dim ) 
```

```{r}
##Integration
features.5000 <- SelectIntegrationFeatures(object.list = Filtered.list, nfeatures = 5000) ##Select features to integrate
Filtered.list <- PrepSCTIntegration(object.list = Filtered.list, anchor.features = features.5000) 
anchors.5000 <- FindIntegrationAnchors(object.list = Filtered.list, normalization.method = "SCT", 
                                       anchor.features = features.5000, dims = 1:10)
combined <- IntegrateData(anchorset = anchors.5000, normalization.method = "SCT", dims = 1:10, new.assay.name = "integrated")
```

```{r}
DefaultAssay(combined) <- "integrated"
set.seed(2022)

##Run dimensionality reduction
combined <- RunPCA(combined)
combined <- RunUMAP(combined, reduction = "pca", dims = 1:10)
combined <- FindNeighbors(combined, reduction = "pca", dims = 1:10)
combined <- FindClusters(combined, resolution = 0.6)
Idents(object = combined) <- "integrated_snn_res.0.6"
DimPlot(combined, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6)
```

```{r}
##Find markers for the cell clusters
##All Markers
All_markers <- FindAllMarkers(object = combined, 
                              only.pos = TRUE,
                              logfc.threshold = 0.25,
                              assay = "RNA")
##Conserved Markers
DefaultAssay(combined) <- "RNA"

C_Markers <- function(cluster){
  FindConservedMarkers(combined,
                       ident.1 = cluster,
                       grouping.var = "Group",
                       only.pos = TRUE) %>%
    rownames_to_column(var = "gene") %>%
    cbind(cluster_id = cluster, .)
}
Conserved_Markers <- map_dfr(c(0:17), C_Markers)

Conserved_Markers <- Conserved_Markers %>% 
  mutate(avg_fc = (Uninfected_avg_log2FC + Infected_avg_log2FC) /2) %>% 
  group_by(cluster_id)
```

```{r}
#Annotation of clusters
cluster.ids <- c("Transitional B cells 2", "Transitional B cells 1", "Transitional B cells 3", "NK cells", "Treg", "T cell 1", "Transitional B cells 4", "CD8a+ T cells", "GC B cells", "Preadipocyte 2", "Plasma cell 1", "Macrophage 1", "T cell 2", "Replicative B cells", "Preadipocyte 1", "Plasma cell 2", "Macrophage 2", "Plasmacytoid DCs")
names(cluster.ids) <- levels(combined)
combined <- RenameIdents(combined, cluster.ids)
DimPlot(combined, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6) + NoLegend() #Figure 6A

#Cell proportions
df <- combined@meta.data %>% group_by(CellType, Group) %>% summarise(Nb = n()) %>% mutate(C = sum(Nb)) %>% mutate(percent = Nb/C*100)
ggplot(df, aes(x = CellType, y = percent, fill = Group))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(percent,"%")), position = position_stack(vjust = 0.5)) #Figure 6B

#Visualise gene expressions
features <- c("Cd79a", "Cd79b", "Ighd", "Jchain", "Cd8a", "Cd3d", "Cd3e", "Trac", "Foxp3", "Icos", "Cd4", "Trbc2", "Ly6d", "Ebf1", "Gzma", "Gzmb", "Nkg7", "Alcda", "Pcna", "Pclaf", "Fn1", "Dcn", "Sparc", "Gsn", "S100a8", "Lyz2", "Fcer1g", "Lgals1", "Ccl5", "S100a4", "Tmem176a", "Tmem176b", "Siglech", "Runx2")
DotPlot(combined, features = features) #Figure 6C
```

##Subsets - T cells
```{r}
DefaultAssay(combined) <- "RNA"
T1 <- subset(combined, idents = c("CD8a+ T cells", "NK cells", "T cell 1", "Treg", "T cell 2")) ##Subset T cell clusters
T1 <- FindVariableFeatures(T1, selection.method = "vst")
T1 <- ScaleData(T1, features = rownames(T1))
DefaultAssay(T1) <- "integrated"

##Run dimensionality reduction
T1 <- RunPCA(T1, features = VariableFeatures(T1), nfeatures.print = 10)
T1 <- RunUMAP(T1, reduction = "pca", dims = 1:10)
T1 <- FindNeighbors(T1, reduction = "pca", dims = 1:10)
T1 <- FindClusters(T1, resolution = 0.3)
Idents(object = T1) <- "integrated_snn_res.0.3"
DimPlot(T1, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6) 

##Find markers
##All Markers
T1_AM <- FindAllMarkers(object = T1, 
                             only.pos = TRUE,
                             logfc.threshold = 0.25,
                             assay = "RNA")  

##Conserved Markers
DefaultAssay(T1) <- "RNA"
T1_CM <- function(cluster){
  FindConservedMarkers(T1,
                       ident.1 = cluster,
                       grouping.var = "Group",
                       only.pos = TRUE) %>%
    rownames_to_column(var = "gene") %>%
    cbind(cluster_id = cluster, .)
}
T1_CM <- map_dfr(c(0:7), T1_CM)
T1_CM <- T1_CM %>% 
  mutate(avg_fc = (Uninfected_avg_log2FC + Infected_avg_log2FC) /2) %>% 
  group_by(cluster_id)

##Annotation of T cell subclusters
T1.ids <- c("Unknown", "Unknown", "Tregs", "CD8+ T cells", "NK cells", "Replicative T cells", "Unknown", "Vgd cells")
names(T1.ids) <- levels(T1)
T1 <- RenameIdents(T1, T1.ids)
DimPlot(T1, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6)

##After removal of unknown clusters
T2 <- subset(T1, idents = c("Tregs", "CD8+ T cells", "NK cells", "Replicative T cells", "Vgd cells")) 
DimPlot(T2, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6, split.by = "Group") #Figure 6D

#Visualise gene expressions
T_features <- c("Icos", "Cd4", "Foxp3", "Nkg7", "Klrd1", "Cd8a", "Cd8b1", "Gzma", "Gzmb", "Top2a", "Mki67", "Hist1h1b", "Cd27", "Tcrg-C1", "Rorc", "Cd163l1", "Il17a")
DotPlot(T2, features = T_features) #Figure 6E
```

```{r}
#IL17ra expression
VlnPlot(combined, features = "Il17ra", split.by = "Group", pt.size = 0) #Figure 7A
```

##Subsets - Adipocytes
```{r}
DefaultAssay(combined) <- "RNA"
A1 <- subset(combined, idents = c("Preadipocyte 1", "Preadipocyte 2")) ##Subset adipocyte clusters
A1 <- FindVariableFeatures(A1, selection.method = "vst")
A1 <- ScaleData(A1, features = rownames(A1))
DefaultAssay(A1) <- "integrated"

##Run dimensionality reduction
A1 <- RunPCA(A1, features = VariableFeatures(A1), nfeatures.print = 10)
A1 <- RunUMAP(A1, reduction = "pca", dims = 1:10)
A1 <- FindNeighbors(A1, reduction = "pca", dims = 1:10)
A1 <- FindClusters(A1, resolution = 0.3)
Idents(object = A1) <- "integrated_snn_res.0.3"
DimPlot(A1, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6) 

##Find markers
##All Markers
A1_AM <- FindAllMarkers(object = A1, 
                             only.pos = TRUE,
                             logfc.threshold = 0.25,
                             assay = "RNA")  

##Conserved Markers
DefaultAssay(A1) <- "RNA"
A1_CM <- function(cluster){
  FindConservedMarkers(A1,
                       ident.1 = cluster,
                       grouping.var = "Group",
                       only.pos = TRUE) %>%
    rownames_to_column(var = "gene") %>%
    cbind(cluster_id = cluster, .)
}
A1_CM <- map_dfr(c(0:4), A1_CM)
A1_CM <- A1_CM %>% 
  mutate(avg_fc = (Uninfected_avg_log2FC + Infected_avg_log2FC) /2) %>% 
  group_by(cluster_id)

##Annotation of T cell subclusters
A1.ids <- c("Interstitial preadipocyte 1", "Committed preadipocyte 1", "Interstitial preadipocyte 2", "Adipogenesis regulatory cells", "Mature adipocytes")
names(A1.ids) <- levels(A1)
A1 <- RenameIdents(A1, A1.ids)
DimPlot(A1, reduction = "umap", label = TRUE, pt.size = 1, label.size = 6, split.by = "Group") #Figure 7B

##Visualise gene expressions
A_features <- c("Ly6a", "Pdgfra", "Cd34", "Dpp4", "Pi16", "Bmp7", "Col4a1", "Col4a2", "Col15a1", "Fabp4", "Plin2", "Cd36", "Fmo2", "F3", "Clec11a", "Pparg")
DotPlot(A1, features = A_features) #Figure 7C

##Cell proportions
df_a1 <- A1@meta.data %>% group_by(Group, CellType) %>% summarise(Nb = n()) %>% mutate(C = sum(Nb)) %>% mutate(proportion = Nb/C)
ggplot(df_a1, aes(x = Group, y = proportion, fill = CellType))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = scales::percent(proportion)), position = position_stack(vjust = 0.5)) #Figure 7D

##IL17ra expression
VlnPlot(A1, features = "Il17ra", split.by = "Group") #Figure 7E

##Lipolysis >20% 
Lipo_high_features <- c("Fabp4", "Pnpla2", "Plaat3", "Mgll")
DotPlot(A1, features = lipo_high_features, split.by = "Group") #Figure 7F

##Lipolysis <20% 
Lipo_low_features <- c("Lipe", "Irs1", "Adrb1", "Ptgs1", "Prkaca", "Plin1")
DotPlot(A1, features = lipo_low_features, split.by = "Group") #Figure 7G
```