neuron_analysis.Rmd

---
title: "analysis"
author: "Yifan Duan"
date: "2024-09-05"
output: html_document
---

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


```{r}
library(Seurat)
library(BPCells)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(patchwork)
library(Matrix)
library(cowplot)
library(irlba)
# set this option when analyzing large datasets
options(future.globals.maxSize = 3e+09)

neuron_obj <- readRDS("data/neuron.Rds")
```


```{r}
# find variable genes with Seurat
neuron_obj <- FindVariableFeatures(neuron_obj, selection.method = "vst", 
                                   nfeatures = 2000, layer = "data")
top20 <- head(VariableFeatures(neuron_obj), 20)

plot1 <- VariableFeaturePlot(neuron_obj)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = TRUE)
plot2

# selecting top 2000 variable genes
neuron_mat <- neuron_obj@assays$RNA$data
neuron_mat_norm <- neuron_mat[VariableFeatures(neuron_obj),]

# saves time according to BPCells
neuron_mat_norm <- neuron_mat_norm %>% write_matrix_dir(tempfile("mat"))
```


```{r}
# PCA
svd <- BPCells::svds(neuron_mat_norm, k = 50)
#svd <- irlba(neuron_mat_norm, nv=50)
pca <- multiply_cols(svd$v, svd$d)

colnames(pca) <- paste0("PC_", 1:ncol(pca))
rownames(pca) <- colnames(neuron_mat_norm)

neuron_obj[["pca"]] <- CreateDimReducObject(embeddings = pca, key = "PC_", 
                                         assay = DefaultAssay(neuron_obj))

plot(svd$d / sqrt(nrow(neuron_mat_norm)))
```


```{r}
# UMAP
set.seed(123)
umap <- uwot::umap(pca)

colnames(umap) <- paste0("UMAP_", 1:2)
rownames(umap) <- colnames(neuron_mat_norm)

neuron_obj[["umap"]] <- CreateDimReducObject(embeddings = umap, key = "UMAP_", 
                                         assay = DefaultAssay(neuron_obj))

DimPlot(neuron_obj, reduction = "umap", group.by = "cell_type__ontology_label")
```


```{r}
# clustering
neuron_obj <- FindNeighbors(neuron_obj, dims = 1:50)
# Define the range of resolutions
resolutions <- c(0.25)

# Loop through the resolutions and store clustering results in Seurat object's metadata
for (res in resolutions) {
  neuron_obj <- FindClusters(neuron_obj, resolution = res, 
                             verbose = FALSE, method = "igraph",
                             algorithm = 4) # leiden
  
  # Store the clustering results with a label corresponding to the resolution
  #cluster_column <- paste0("RNA_snn_res.", res)
}

DimPlot(neuron_obj, reduction = "umap", group.by = "RNA_snn_res.0.25", label = T) + 
  theme(legend.position = "none")
DimPlot(neuron_obj, reduction = "umap", group.by = "cell_type_custom", label = T) + 
  theme(legend.position = "none")

```


```{r}
Idents(neuron_obj) <- "RNA_snn_res.0.25"
#FindMarkers(neuron_obj, ident.1 = "7", assay = "RNA")

dim(neuron_obj[["RNA"]]$counts)

pdf("figure/marker_mature_OSN.pdf")
FeaturePlot(neuron_obj, features = c("Omp", "Stoml3", "Cnga2", "Adcy3"))
dev.off()
```


```{r investigating doublet}
# 1. higher counts than the rest
# 2. mixture of markers 
grp7_marker <- FindMarkers(neuron_obj, ident.1 = "7", assay = "RNA")

library(EnhancedVolcano)
pdf("figure/DE_genes_grp7.pdf")
EnhancedVolcano(grp7_marker, 
                rownames(grp7_marker),
                x ="avg_log2FC", 
                y ="p_val_adj")
dev.off()


VlnPlot(neuron_obj, features = "number_of_reads", group.by = "seurat_clusters")
VlnPlot(neuron_obj, features = "number_of_features", group.by = "seurat_clusters")
```


```{r}
saveRDS(neuron_obj, file = "data/neuron.Rds")

FeaturePlot(neuron_obj, features = c("Nqo1","Acsm4", "Nfix", "Ncam2"))
# dorsal acsm4, nqo1
# ventral nfix, ncam2
```