---
title: "Subclustering astro"
author: "Shahzad"
date: '2023-11-01'
output: html_document
---

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

```{r}
#Create astrocyte object (astro) for subcluster analysis

astro <- (subset(gswt, idents =c("Astrocyte")))
astro
```

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

```{r}
#Specify feature minimum and maximum, and exclude high percent.mt (above 2.5)
astro <- subset(astro, subset = nFeature_RNA > 200 & nFeature_RNA < 2000 & percent.mt < 2.5 )

#Re-visualize astro object features, counts, and percent.mt
VlnPlot(astro, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

```

```{r}
#Perform normalization, find variable features, scale, run PCA, and find subclusters
astro <- NormalizeData(astro, normalization.method = "LogNormalize", scale.factor = 10000)

astro <- FindVariableFeatures(astro, selection.method = "vst", nfeatures = 2500)

astro_top10 <- head(VariableFeatures(astro), 10)

plot1 <- VariableFeaturePlot(astro)
plot2 <- LabelPoints(plot = plot1, points = astro_top10, repel = TRUE)
plot1 + plot2

all.genes <- rownames(astro)
astro <- ScaleData(astro)

astro <- RunPCA(astro, features = VariableFeatures(object = astro))

astro <- FindNeighbors(astro, dims = 1:50)
astro <- FindClusters(astro, resolution = 0.5)
```

```{r}
#Run TSNE and visualize TSNE astro_plot
astro <- RunTSNE(astro, dims = 1:50)
astro_plot <- DimPlot(astro, reduction = "tsne", label = TRUE)
astro_plot
```

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

```{r}
#Remove Fibrolast-like subcluster 5
astro <- subset(astro, idents =c("5"), invert = TRUE)
```

```{r}
#Normalize, Find variable features, scale data, run PCA, and find subclusters 
astro <- NormalizeData(astro, normalization.method = "LogNormalize", scale.factor = 10000)

astro <- FindVariableFeatures(astro, selection.method = "vst", nfeatures = 2500)

astro_top10 <- head(VariableFeatures(astro), 10)

plot1 <- VariableFeaturePlot(astro)
plot2 <- LabelPoints(plot = plot1, points = astro_top10, repel = TRUE)
plot1 + plot2

all.genes <- rownames(astro)
astro <- ScaleData(astro)

astro <- RunPCA(astro, features = VariableFeatures(object = astro))

astro <- FindNeighbors(astro, dims = 1:50)
astro <- FindClusters(astro, resolution = 0.5)
```

```{r}
#Run TSNE and visualize TSNE astro_plot
astro <- RunTSNE(astro, dims = 1:50)
astro_plot <- DimPlot(astro, reduction = "tsne", label = TRUE)
astro_plot
```

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

```{r}
#Remove non-astrocytic subcluster 5
astro <- subset(astro, idents =c("5"), invert = TRUE)
```

```{r}
#Normalize, Find variable features, scale data, run PCA, and find subclusters 
astro <- NormalizeData(astro, normalization.method = "LogNormalize", scale.factor = 10000)

astro <- FindVariableFeatures(astro, selection.method = "vst", nfeatures = 2500)

astro_top10 <- head(VariableFeatures(astro), 10)

plot1 <- VariableFeaturePlot(astro)
plot2 <- LabelPoints(plot = plot1, points = astro_top10, repel = TRUE)
plot1 + plot2

all.genes <- rownames(astro)
astro <- ScaleData(astro)

astro <- RunPCA(astro, features = VariableFeatures(object = astro))

astro <- FindNeighbors(astro, dims = 1:50)
astro <- FindClusters(astro, resolution = 0.5)
```

```{r}
#Run TSNE and visualize TSNE astro_plot
astro <- RunTSNE(astro, dims = 1:50)
astro_plot <- DimPlot(astro, reduction = "tsne", label = TRUE)
astro_plot

```

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

```{r}
#Assign new cluster IDs based on Top 10 markers
new.cluster.ids <- c("Aldh1a1", "Crym", "Crym", "Nrg1", "Gfap")
names(new.cluster.ids) <- levels(astro)
astro <- RenameIdents(astro, new.cluster.ids)

#Visualize astro plot with new cluster IDs
astro_plot <- DimPlot(astro, reduction = "tsne", label = TRUE , pt.size = 0.5) + NoLegend()
astro_plot
```

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

```{r}
#Create columns containing genotype information
astro <-AddMetaData(astro, metadata = substr(astro@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 = gswt, ident = "INSERT_CLUSTER_ID_NAME"),names(which(gswt$genotype=='GS')))
wt_cellname <-intersect(WhichCells(object = gswt, ident = "INSERT_CLUSTER_ID_NAME"),names(which(gswt$genotype=='WT')))
```

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