---
title: "Chat Analysis"
author: "Shahzad Khan"
date: '2023-11-01'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library(dplyr)
library(Seurat)
library(patchwork)
```

```{r}
#Create chat object for subcluster anaylsis
chat <- (subset(gswt, idents =c("Chat_IN")))
chat
```

```{r}
#Visualize features, counts, and percent.mt
VlnPlot(chat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), slot = "counts", ncol = 3)
```

```{r}
#Normalize
chat <- NormalizeData(chat, normalization.method = "LogNormalize", scale.factor = 10000)
```

```{r}
#Find variable features
chat <- FindVariableFeatures(chat, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(chat), 10)

#Plot variable features with and without labels
plot1 <- VariableFeaturePlot(chat)
plot2 <- LabelPoints(plot = plot1, points = top10)
plot1 + plot2
```

```{r}
#Scale object
all.genes <- rownames(chat)
chat <- ScaleData(chat, features = all.genes)
```

```{r}
#Perform PCA
chat <- RunPCA(chat, features = VariableFeatures(object = chat))
```

```{r}
#Find subclusters
chat <- FindNeighbors(chat, dims = 1:50)
chat <- FindClusters(chat, resolution = .5)
```

```{r}
#Run TSNE and visualize chat_plot
chat <- RunTSNE(chat, dims = 1:50)
chat_plot <- DimPlot(chat, reduction = "tsne", label.size = 5, label = TRUE)
chat_plot
```

```{r}
#Find markers for chat subclusters
chat.markers <- FindAllMarkers(chat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
chat.markers %>% 
    group_by(cluster) %>%
    slice_max(n = 10, order_by = avg_log2FC)
```

```{r}
#Remove oligodendrocyte-like subcluster 2
chat <- subset(chat, idents =c("2"), invert = TRUE)
chat
```

```{r}
#Normalize again after removing subcluster 2
chat <- NormalizeData(chat, normalization.method = "LogNormalize", scale.factor = 10000)
```

```{r}
#Find variable features
chat <- FindVariableFeatures(chat, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(chat), 10)

#Plot variable features with and without labels
plot1 <- VariableFeaturePlot(chat)
plot2 <- LabelPoints(plot = plot1, points = top10)
plot1 + plot2
```

```{r}
#Scale object
all.genes <- rownames(chat)
chat <- ScaleData(chat, features = all.genes)
```

```{r}
#Perform PCA
chat <- RunPCA(chat, features = VariableFeatures(object = chat))
```

```{r}
#Find subclusters
chat <- FindNeighbors(chat, dims = 1:50)
chat <- FindClusters(chat, resolution = .5)
```

```{r}
#Run TSNE and visualize chat_plot
chat <- RunTSNE(chat, dims = 1:50)
chat_plot <- DimPlot(chat, reduction = "tsne", label.size = 5, label = TRUE)
chat_plot
```

```{r}
#Find markers for chat subclusters
chat.markers <- FindAllMarkers(chat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
chat.markers %>% 
    group_by(cluster) %>%
    slice_max(n = 10, order_by = avg_log2FC)

```

```{r}
#Rename subclusters and visualize
new.cluster.ids <- c("Elavl2","Grm5")
names(new.cluster.ids) <- levels(chat)
chat <- RenameIdents(chat, new.cluster.ids)
DimPlot(chat, reduction = "tsne", label = TRUE, pt.size = 0.5) + NoLegend()
```

```{r}
#Find markers for chat subclusters using new cluster IDs
chat.markers <- FindAllMarkers(chat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
chat.markers %>% 
    group_by(cluster) %>%
    slice_max(n = 10, order_by = avg_log2FC)
```


```{r}
#Visualize gene-specific expression using feature plot
FeaturePlot(chat, features = "GENE_OF_INTEREST", reduction = "tsne", label = TRUE)
```

```{r}
#Create columns containing genotype information
chat <-AddMetaData(chat, metadata = substr(chat@meta.data$orig.ident, 1,2), col.name = "genotype")
```

```{r}
#Isolate information for G2019S nuclei (gs_cellname) and WT nuclei (wt_cellname)
gs_cellname <-intersect(WhichCells(object = chat, ident = "INSERT_CLUSTER_ID_NAME"),names(which(chat$genotype=='GS')))
wt_cellname <-intersect(WhichCells(object = chat, ident = "INSERT_CLUSTER_ID_NAME"),names(which(chat$genotype=='WT')))
```

```{r}
#Compare gene expression features of G2019S nuclei (gs_cellname) to WT nuclei (wtcellname) and visualize
comp_markers <- FindMarkers(chat, ident.1 = gs_cellname, ident.2 = wt_cellname, logfc.threshold = 0)
comp_markers
```

