---
title: "10x snRNAseq pipeline"
author: "Shahzad Khan"
date: "2023-11-01"
output: html_document
---

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.

```{r}
library(dplyr)
library(Seurat)
library(patchwork)
library(ggplot2)
install.packages("svglite")
suppressMessages(require(Seurat))
suppressMessages(require(Matrix))
install.packages("remotes")
remotes::install_github("chris-mcginnis-ucsf/DoubletFinder", upgrade = F, force = TRUE)
suppressMessages(require(DoubletFinder))
```

```{r}
#Load the GSA dataset
  gsa.data <- Read10X(data.dir = "10X/GSA/")
gsa <- CreateSeuratObject(counts = gsa.data, project = "GSA", min.cells = 3, min.features = 200)
gsa
```

```{r}
#Add percent feature
gsa[["percent.mt"]] <- PercentageFeatureSet(gsa, pattern = "^mt-")
```

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

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

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

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

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

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

```{r}
#Find Clusters
gsa <- FindNeighbors(gsa, dims = 1:40)
gsa <- FindClusters(gsa, resolution = 0.5)
```

```{r}
#RunUMAP
gsa <- RunUMAP(gsa, dims = 1:30)
```

```{r}
#Visualize clusters
DimPlot(gsa, reduction = "umap", label.size = 3, label = TRUE)
```

```{r}
#Identify high "percent.mt" cluster(s)
FeaturePlot(gsa, features = "percent.mt", label=T)
```

```{r}
#Remove high "percent.mt" cluster(s)
gsa <- subset(gsa, idents =c("8"), invert = TRUE)
```

```{r}
#Identify cluster(s) for doublet finder 
gsa = FindVariableFeatures(gsa, verbose = F)
gsa = ScaleData(gsa, vars.to.regress = c("nFeature_RNA", "percent.mt"),
    verbose = F)
gsa = RunPCA(gsa, verbose = F, npcs = 50)
gsa = RunUMAP(gsa, dims = 1:50, verbose = F)
```

```{r}
#Run parameter optimization with paramSweep
sweep.res <- paramSweep_v3(gsa) 
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) 
bcmvn <- find.pK(sweep.stats)
barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2, axes = TRUE)
```

```{r}
nExp <- round(ncol(gsa) * 0.039)  # expect 3.9% doublets
gsa <- doubletFinder_v3(gsa, pN = 0.25, pK = 0.005, nExp = nExp, PCs = 1:10)
```

```{r}

#Name of the DF prediction can change, so extract the correct column name.
DF.name = colnames(gsa@meta.data)[grepl("DF.classification", colnames(gsa@meta.data))]

cowplot::plot_grid(ncol = 2, DimPlot(gsa, group.by = "orig.ident") + NoAxes(),
    DimPlot(gsa, group.by = DF.name) + NoAxes())
```

```{r}
#Remove predicted doublets; keep singlets
gsa = gsa[, gsa@meta.data[, DF.name] == "Singlet"]
dim(gsa)
```

```{r}
#Load the GSB dataset
  gsb.data <- Read10X(data.dir = "10X/GSB/")
gsb <- CreateSeuratObject(counts = gsb.data, project = "GSB", min.cells = 3, min.features = 200)
gsb
```

```{r}
#Add percent feature
gsb[["percent.mt"]] <- PercentageFeatureSet(gsb, pattern = "^mt-")
```

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

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

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

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

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

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

```{r}
#Find Clusters
gsb <- FindNeighbors(gsb, dims = 1:40)
gsb <- FindClusters(gsb, resolution = 0.5)
```

```{r}
#RunUMAP
gsb <- RunUMAP(gsb, dims = 1:30)
```

```{r}
#Visualize clusters
DimPlot(gsb, reduction = "umap", label.size = 3, label = TRUE)
```

```{r}
#Identify high "percent.mt" cluster(s)
FeaturePlot(gsb, features = "percent.mt", label=T)
```

```{r}
#Remove high "percent.mt" cluster(s)
gsb <- subset(gsb, idents =c("14"), invert = TRUE)
```

```{r}
#Identify cluster(s) for doublet finder
gsb = FindVariableFeatures(gsb, verbose = F)
gsb = ScaleData(gsb, vars.to.regress = c("nFeature_RNA", "percent.mt"),
    verbose = F)
gsb = RunPCA(gsb, verbose = F, npcs = 50)
gsb = RunUMAP(gsb, dims = 1:50, verbose = F)
```

```{r}
# Run parameter optimization with paramSweep
sweep.res <- paramSweep_v3(gsb) 
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) 
bcmvn <- find.pK(sweep.stats)
barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2, axes = TRUE)
```

```{r}
nExp <- round(ncol(gsb) * 0.05)  # expect 5% doublets
gsb <- doubletFinder_v3(gsb, pN = 0.25, pK = 0.005, nExp = nExp, PCs = 1:10)
```

```{r}
#Name of the DF prediction can change, so extract the correct column name.
DF.name = colnames(gsb@meta.data)[grepl("DF.classification", colnames(gsb@meta.data))]



cowplot::plot_grid(ncol = 2, DimPlot(gsb, group.by = "orig.ident") + NoAxes(),
    DimPlot(gsb, group.by = DF.name) + NoAxes())
```

```{r}
#Remove predicted doublets; keep singlets
gsb = gsb[, gsb@meta.data[, DF.name] == "Singlet"]
dim(gsb)
```

```{r}
#Load the GSC dataset
  gsc.data <- Read10X(data.dir = "10X/GSC/")
gsc <- CreateSeuratObject(counts = gsc.data, project = "GSC", min.cells = 3, min.features = 200)
gsc
```

```{r}
#Add percent feature
gsc[["percent.mt"]] <- PercentageFeatureSet(gsc, pattern = "^mt-")
```

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

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

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

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

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

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

```{r}
#Find Clusters
gsc <- FindNeighbors(gsc, dims = 1:40)
gsc <- FindClusters(gsc, resolution = 0.5)
```

```{r}
#RunUMAP
gsc <- RunUMAP(gsc, dims = 1:30)
```

```{r}
#Visualize clusters
DimPlot(gsc, reduction = "umap", label.size = 3, label = TRUE)
```

```{r}
#Identify high "percent.mt" cluster(s)
FeaturePlot(gsc, features = "percent.mt", label=T)
```

```{r}
#Remove high "percent.mt" cluster(s)
gsc <- subset(gsc, idents =c("3","16"), invert = TRUE)
```

```{r}
#Identify cluster(s) for doublet finder
gsc = FindVariableFeatures(gsc, verbose = F)
gsc = ScaleData(gsc, vars.to.regress = c("nFeature_RNA", "percent.mt"),
    verbose = F)
gsc = RunPCA(gsc, verbose = F, npcs = 50)
gsc = RunUMAP(gsc, dims = 1:50, verbose = F)
```

```{r}
#Run parameter optimization with paramSweep
sweep.res <- paramSweep_v3(gsc) 
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) 
bcmvn <- find.pK(sweep.stats)
barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2, axes = TRUE)
```

```{r}
nExp <- round(ncol(gsc) * 0.05)  # expect 4% doublets
gsc <- doubletFinder_v3(gsc, pN = 0.25, pK = 0.001, nExp = nExp, PCs = 1:10)
```

```{r}
#Name of the DF prediction can change, so extract the correct column name.
DF.name = colnames(gsc@meta.data)[grepl("DF.classification", colnames(gsc@meta.data))]



cowplot::plot_grid(ncol = 2, DimPlot(gsc, group.by = "orig.ident") + NoAxes(),
    DimPlot(gsc, group.by = DF.name) + NoAxes())
```

```{r}
#Remove predicted doublets; keep singlets
gsc = gsc[, gsc@meta.data[, DF.name] == "Singlet"]
dim(gsc)
```


```{r}
#Load the WTA dataset
  wta.data <- Read10X(data.dir = "10X/WTA/")
wta <- CreateSeuratObject(counts = wta.data, project = "WTA", min.cells = 3, min.features = 200)
wta
```

```{r}
#Add percent feature
wta[["percent.mt"]] <- PercentageFeatureSet(wta, pattern = "^mt-")
```

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

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

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

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

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

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

```{r}
#Find Clusters
wta <- FindNeighbors(wta, dims = 1:40)
wta <- FindClusters(wta, resolution = 0.5)
```

```{r}
#RunUMAP
wta <- RunUMAP(wta, dims = 1:30)
```

```{r}
#Visualize clusters
DimPlot(wta, reduction = "umap", label.size = 3, label = TRUE)
```

```{r}
#Identify high "percent.mt" cluster(s)
FeaturePlot(wta, features = "percent.mt", label=T)
```

```{r}
#Remove high "percent.mt" cluster(s)
wta <- subset(wta, idents =c("9"), invert = TRUE)
```

```{r}
#Identify cluster(s) for doublet finder
wta = FindVariableFeatures(wta, verbose = F)
wta = ScaleData(wta, vars.to.regress = c("nFeature_RNA", "percent.mt"),
    verbose = F)
wta = RunPCA(wta, verbose = F, npcs = 50)
wta = RunUMAP(wta, dims = 1:50, verbose = F)
```

```{r}
#Run parameter optimization with paramSweep
sweep.res <- paramSweep_v3(wta) 
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) 
bcmvn <- find.pK(sweep.stats)
barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2, axes = TRUE)
```


```{r}
nExp <- round(ncol(wta) * 0.046)  # expect 4.6% doublets
wta <- doubletFinder_v3(wta, pN = 0.25, pK = 0.29, nExp = nExp, PCs = 1:10)
```

```{r}
#Name of the DF prediction can change, so extract the correct column name.
DF.name = colnames(wta@meta.data)[grepl("DF.classification", colnames(wta@meta.data))]



cowplot::plot_grid(ncol = 2, DimPlot(wta, group.by = "orig.ident") + NoAxes(),
    DimPlot(wta, group.by = DF.name) + NoAxes())
```

```{r}
#Remove predicted doublets; keep singlets
wta = wta[, wta@meta.data[, DF.name] == "Singlet"]
dim(wta)
```

```{r}
#Load the WTB dataset
  wtb.data <- Read10X(data.dir = "10X/WTB/")
wtb <- CreateSeuratObject(counts = wtb.data, project = "WTB", min.cells = 3, min.features = 200)
wtb
```

```{r}
#Add percent feature
wtb[["percent.mt"]] <- PercentageFeatureSet(wtb, pattern = "^mt-")
```

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

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

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

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

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

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

```{r}
#Find Clusters
wtb <- FindNeighbors(wtb, dims = 1:40)
wtb <- FindClusters(wtb, resolution = 0.5)
```

```{r}
#RunUMAP
wtb <- RunUMAP(wtb, dims = 1:30)
```

```{r}
#Visualize clusters
DimPlot(wtb, reduction = "umap", label.size = 3, label = TRUE)
```

```{r}
#Identify high "percent.mt" cluster(s)
FeaturePlot(wtb, features = "percent.mt", label=T)
```

```{r}
#Remove high "percent.mt" cluster(s)
wtb <- subset(wtb, idents =c("16"), invert = TRUE)
```

```{r}
#Identify cluster(s) for doublet finder
wtb = FindVariableFeatures(wtb, verbose = F)
wtb = ScaleData(wtb, vars.to.regress = c("nFeature_RNA", "percent.mt"),
    verbose = F)
wtb = RunPCA(wtb, verbose = F, npcs = 50)
wtb = RunUMAP(wtb, dims = 1:50, verbose = F)
```

```{r}
#Run parameter optimization with paramSweep
sweep.res <- paramSweep_v3(wtb) 
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) 
bcmvn <- find.pK(sweep.stats)
barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2, axes = TRUE)
```

```{r}
nExp <- round(ncol(wtb) * 0.05)  # expect 5% doublets
wtb <- doubletFinder_v3(wtb, pN = 0.25, pK = 0.01, nExp = nExp, PCs = 1:10)
```

```{r}
#Name of the DF prediction can change, so extract the correct column name.
DF.name = colnames(wtb@meta.data)[grepl("DF.classification", colnames(wtb@meta.data))]



cowplot::plot_grid(ncol = 2, DimPlot(wtb, group.by = "orig.ident") + NoAxes(),
    DimPlot(wtb, group.by = DF.name) + NoAxes())
```

```{r}
#Remove predicted doublets; keep singlets
wtb = wtb[, wtb@meta.data[, DF.name] == "Singlet"]
dim(wtb)
```

```{r}
#Load the WTC dataset
  wtc.data <- Read10X(data.dir = "10X/WTC/")
wtc <- CreateSeuratObject(counts = wtc.data, project = "WTC", min.cells = 3, min.features = 200)
wtc
```

```{r}
#Add percent feature
wtc[["percent.mt"]] <- PercentageFeatureSet(wtc, pattern = "^mt-")
```

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

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

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

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

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

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

```{r}
#Find Clusters
wtc <- FindNeighbors(wtc, dims = 1:40)
wtc <- FindClusters(wtc, resolution = 0.5)
```

```{r}
#RunUMAP
wtc <- RunUMAP(wtc, dims = 1:30)
```

```{r}
#Visualize clusters
DimPlot(wtc, reduction = "umap", label.size = 3, label = TRUE)
```

```{r}
#Identify high "percent.mt" cluster(s)
FeaturePlot(wtc, features = "percent.mt", label=T)
```

```{r}
#Remove high "percent.mt" cluster(s)
wtc <- subset(wtc, idents =c("13"), invert = TRUE)
```

```{r}
#Identify cluster(s) for doublet finder
wtc = FindVariableFeatures(wtc, verbose = F)
wtc = ScaleData(wtc, vars.to.regress = c("nFeature_RNA", "percent.mt"),
    verbose = F)
wtc = RunPCA(wtc, verbose = F, npcs = 50)
wtc = RunUMAP(wtc, dims = 1:50, verbose = F)
```

```{r}
# Run parameter optimization with paramSweep
sweep.res <- paramSweep_v3(wtc) 
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) 
bcmvn <- find.pK(sweep.stats)
barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2, axes = TRUE)
```


```{r}
nExp <- round(ncol(wtc) * 0.061)  # expect 6.1% doublets
wtc <- doubletFinder_v3(wtc, pN = 0.25, pK = 0.04, nExp = nExp, PCs = 1:10)
```

```{r}
#Name of the DF prediction can change, so extract the correct column name.
DF.name = colnames(wtc@meta.data)[grepl("DF.classification", colnames(wtc@meta.data))]



cowplot::plot_grid(ncol = 2, DimPlot(wtc, group.by = "orig.ident") + NoAxes(),
    DimPlot(wtc, group.by = DF.name) + NoAxes())
```

```{r}
#Remove predicted doublets; keep singlets
wtc = wtc[, wtc@meta.data[, DF.name] == "Singlet"]
dim(wtc)
```

```{r}
#Merge gsa and gsb (creates gs1)
gs1 <- merge(gsa, y = gsb, add.cell.ids = c("gsa", "gsb"), project = "gs1")
gs1
```

```{r}
#Merge gs1 and gsc (creates gs2)
gs2 <- merge(gs1, y = gsc, add.cell.ids = c("gs1", "gsc"), project = "gs2")
gs2
```

```{r}
#Merge wta and wtb (creates wt1)
wt1 <- merge(wta, y = wtb, add.cell.ids = c("wta", "wtb"), project = "wt1")
wt1
```

```{r}
#Merge wt1 and wtc (creates wt2)
wt2 <- merge(wt1, y = wtc, add.cell.ids = c("wt1", "wtc"), project = "wt2")
wt2
```

```{r}
#Merge gs2 and wt2 (creates gswt)
gswt <- merge(gs2, y = wt2, add.cell.ids = c("gs2", "wt2"), project = "gswt")
gswt
```

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

```{r}
#Add percent feature
gswt <- FindVariableFeatures(gswt, selection.method = "vst", nfeatures = 2500)

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

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

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

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


```{r}
#Find Clusters
gswt <- FindNeighbors(gswt, dims = 1:40)
gswt <- FindClusters(gswt, resolution = 0.5)
```

```{r}
#Run UMAP
gswt <- RunUMAP(gswt, dims = 1:40)
```

```{r}
#Visualize clusters
DimPlot(gswt, reduction = "umap", label.size = 3, label = TRUE)
```

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

```{r}
#Remove high percent mt- cluster 27 using Top10 marker output
gswt <- subset(gswt, idents =c("27"), invert = TRUE)
```

```{r}
#Normalize gswt after cluster 27 removal
gswt <- NormalizeData(gswt, normalization.method = "LogNormalize", scale.factor = 10000)
```

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

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

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

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

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

```{r}
#Find clusters
gswt <- FindNeighbors(gswt, dims = 1:40)
gswt <- FindClusters(gswt, resolution = 0.5)
```

```{r}  
#Run UMAP
gswt <- RunUMAP(gswt, dims = 1:40)
```

```{r}
#Visualize Clusters
DimPlot(gswt, reduction = "umap", label.size = 3, label = TRUE)
```

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

```{r}
#Rename clusters based on Top 10 Markers and known cell specific markers
new.cluster.ids <- c("dSPN", "iSPN", "Oligo", "Astrocyte", "dSPN", "iSPN", "Microglia", "Oligo", "NPC", "eSPN", "Astrocyte", "Pvalb_IN", "Endothelial", "Sst_IN", "Ependyma", "Mural", "Chat_IN", "Microglia", "OPC")
names(new.cluster.ids) <- levels(gswt)
gswt <- RenameIdents(gswt, new.cluster.ids)
```

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

```{r}
#Run TSNE and visualize clusters
gswt <- RunTSNE(gswt, dims = 1:40)
DimPlot(gswt, reduction = "tsne", label = TRUE, pt.size = 0.5)
```

```{r}
#Create columns containing genotype information
gswt <-AddMetaData(gswt, metadata = substr(gswt@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
#This and the preceding command are used to compare genotype-specific differences between G2019S and WT nuclei
comp_markers <- FindMarkers(gswt, ident.1 = gs_cellname, ident.2 = wt_cellname, logfc.threshold = 0)
comp_markers
```

```{r}
#Save gene markers for gswt
saveRDS(gswt.markers, file = "gswt.markers_with_cluster_IDs")
```

```{r}
#Save gswt object with cluster IDs
gswt <- readRDS(file ="gswt_with_cluster_IDs")
```

```{r}
#Create Astrocyte object (astro) for subcluster analysis
astro <- subset(gswt, idents =c("Astrocyte"))
astro
```

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