# Set up

To analyze single nuclei/cell data you will most likely need specialized
software (aka. "packages"). We have already installed what you need for
this tutorial, so you just need to tell R where to find it and "load"
it. The `.libPaths` function helps you see or change the locations where R
will look for software, we will use it to include a new location
stored in the `tutorial_packages` variable, then we will load the packages we 
will use:

```{r libs, warning=FALSE, message=FALSE}
tutorial_packages <- "/common/tutorials/bioinformatics/single_cell/software/r-rich4.6_lib/" 
.libPaths(c(tutorial_packages, .libPaths()))
.libPaths()

library(ggplot2)
library(tidyverse)
library(Seurat)

```

```{r chunck_opts, include=FALSE, warning=FALSE, error=FALSE}

knitr::opts_chunk$set(warning = FALSE, message = FALSE, fig.align='center')

```

When you generate results, they will be stored by default in your
current working directory, so it is handy to learn how to find or change
it. This uses the `getwd` and `setwd` functions.

```{r setup, warning=FALSE}

# Here is where your outputs will be stored by default:
getwd()

# If you want to change it:
project_dir <- "/projects/standard/riss/calix015/sc_tutorial/"
setwd(project_dir)

```

Finally, for this tutorial we will be using data downloaded from GEO
with the accession GSE247974. Note that we are assigning the location of
this data (aka. its "path") to the variable `data_dir` so that we can
reuse it easily.

```{r data}

data_dir <- "/common/tutorials/bioinformatics/single_cell/data/neuro_dataset/" 

# You can see the contents of this directory with the "list.files" function:
list.files(data_dir)
list.files(paste0(data_dir,"GSE247974_D3"))

```

# How to load the count data?

This depends on the format of the input. The `Seurat` package has
different functions for different types of data and you can browse and find 
more information [here](https://satijalab.org/seurat/reference/). For this 
data set (10X) we will use the `Read10X` function to store the data in the 
`so` variable:

```{r load_mtx}

count_mtx <- Read10X(data.dir = paste0(data_dir,"GSE247974_D3"))

```

We can see details about how to use this function (or any other) and the default 
parameters with the `?` command:

```{r help, eval=FALSE}

?Read10X

```

Now we will inspect the size and contents of the count matrix we just loaded. 
We will use a subset of the rows and columns to get a little preview instead 
of showing the entire matrix with `[ , ]`:

```{r sparse_mtx}

dim(count_mtx)
count_mtx[1:10,1:5]

```

This small example we have chosen for the tutorial has information from 36,601 
genes and 4,245 barcodes. Note that most of the values are `.`, which represent 
zeros. This is what a "sparse" data set looks like. Also note that the row names 
are gene names and the column names are barcodes.


# Preprocessing overview

## Format the counts matrix

To analyze a data set like this we should first format it so that the packages 
we will use can access the information. The `CreateSeuratObject` function does that:

```{r fmt}
so <- CreateSeuratObject(count_mtx)
so
```
Since we know there are a lot of zeros, we can use some minimum thresholds
to filter out the barcodes that are most likely empty droplets. 
We will rewrite the `so` variable:

```{r so_filter}

# Create a SO with filters
so <- CreateSeuratObject(counts = count_mtx, min.features = 10, min.cells = 5)
so

```

## Structure of a Seurat object (so)

You can browse the structure of the data by clicking on its name in the `Data` 
panel or with the `View` function. The `head` function shows you a little preview 
instead of the entirety of the data.
We will peek into a few key elements:

```{r so_struct}
# all the genes 
rownames(so) %>% head()

# the counts
so@assays$RNA$counts[1:10,1:5]

# the metadata
so@meta.data %>% head()

```

## QC to clustering

Please check out the Seurat [tutorials](https://satijalab.org/seurat/articles/pbmc3k_tutorial)
for detailed information. Briefly, this process involves a quality control check 
and filtering out barcodes that don't contain acceptable cells/nuclei, 
count normalization, exploring the variability of the data, scaling, 
dimensionality reduction, clustering, cluster identification and sometimes sample
integration.

A template workflow would look like this:

```{r preproc, eval=FALSE}

# calculate filtering metrics:
so[["percent.mt"]] <- PercentageFeatureSet(so, pattern = "^MT-")
# plot values and select thresholds:
VlnPlot(so,c("nCount_RNA", "nFeature_RNA","percent.mt"))
# filter:
so <- subset(so, 
             subset = nFeature_RNA > 1000 & nFeature_RNA < 7500 & nCount_RNA < 50000)
# Normalize:
so <- NormalizeData(so) 
# Inspect variability:
so <- FindVariableFeatures(so) 
# Scale:
so <- ScaleData(so, features = rownames(so)) 
# Reduce dimensionality:
so <- RunPCA(so, 
              features = VariableFeatures(object = so)) 
# Define dimensionality threshold:
ElbowPlot(so) 
ndims <- 15
# Cluster:
so <- FindNeighbors(so, dims = 1:ndims) 
so <- FindClusters(so, resolution = 0.5) 
# Calculate other data projections for plotting:
so <- RunUMAP(so, dims = 1:ndims)

```

But be aware that the thresholds and options chosen in these functions
should mold to your data's characteristics. 
*Please think carefully about these steps because they are the base for the results you will get.*

Finally, for cluster or cell identification we need information about the 
samples and the expected cell types, marker genes or a really good database
to compare the data against. This process is quite particular to your samples,
iterative and sometimes quite difficult. 
We talk a little more in depth about it in our "Introduction to Single Cell 
Genomics" tutorial.

# Pre-analyzed object 

In this part we will show you how to open a pre-analyzed object in case you need 
to, for example, generate plots of your favorite genes or cell types.

```{r pre_analized}
# load the pre-processed and labeled data
so_done <- readRDS(paste0(data_dir,"GSE247974_D3.Rds"))
so_done
```

After analysis we have more layers with gene expression information, a set of 
variable genes, dimensional reductions and more information in the metadata:

```{r meta}
so_done@meta.data %>%head()
```

The object meta data is nothing more than a data frame with the cell barcodes as 
row names and different cell attributes in the columns, so we can count and 
summarize this information as needed. For example, we can look at the "celltype"
column and count how many cells of each type are present:

```{r cell_n}
table(so_done$celltype) 
```

To assign each cell an identity value stored in a metadata column, we use the 
`Idents` function:

```{r ident}
Idents(so_done) <- so_done$celltype
```


## Subset cells

You can also use the metadata to subset a group of cells that most interest you. 
For example, let's subset all cell types that have at least 10 representatives:

```{r subset}
# select the ids that you want with a condition:
ids_ok <- table(so_done$celltype) > 10
ids_ok <- names(ids_ok)[ids_ok]
#subset:
so_done_subset <- subset(so_done, idents = ids_ok)
```


## Plot cell projections

The fisrt 2 dimensions after the dimensionality reduction step are usually 
chosen for a 2D representation of your data. Commonly, a PCA reduction is the 
basis for all downstream analysis, let's see how the cells look when projected 
in a 2D space:

```{r pca}
DimPlot(so_done_subset, reduction = "pca")
```

Often, you'll see another very common projection strategy called UMAP:

```{r umap}
DimPlot(so_done_subset)
```
Note that neither of these plots tells you the full story of the data, it's 
just a "summary picture" and a very common mistake is to misinterpret the 
distance between clusters in a UMAP with the difference between cells. 
This is generally not correct and can lead you to wrong conclusions. 

You can customize your plot changing the function parameters (use `?DimPlot`). 
Here I am removing the legend and annotating the cell type over the plot:

```{r umap_custom}
DimPlot(so_done_subset,label = TRUE, pt.size = 0.5) + NoLegend()
```


## Plot gene expression

Two very common plotting options to summarize the gene expression of your genes of 
interest across different samples or cell groups are the violin and dot plots. 
Let's make some with the most variable genes in this data set:

```{r gene}
# top 5 var genes
top5 <- so_done_subset@assays$SCT@var.features[1:5]

# Violin of the most variable:
VlnPlot(so_done_subset, top5[1])+NoLegend()+xlab("")

# Dotplot:
DotPlot(so_done_subset, top5)
```

Let's improve this dotplot a little:

```{r custom_dot}
holder <- DotPlot(so_done_subset, top5)
holder <- holder+ggtitle("customizing my plot")+
  scale_color_gradient(low = "grey",high = "red")+
  xlab("")+
  ylab("")+
  theme(axis.text = element_text(size=8))
holder

```

If you'd like to see a particular gene's expression data or a metadata column
with quantitative information projected on the UMAP plot, you can use a 
"feature" plot:

```{r feat}
colnames(so_done_subset@meta.data)

# Gene expression:
FeaturePlot(so_done_subset,  top5[1])

# Metadata column:
FeaturePlot(so_done_subset, "nCount_SCT")
```

# Save your work!

Remember to save your results and plots before your session ends. 

```{r save}
# Save the SO:
saveRDS(so_done_subset, paste0(project_dir,"subset.rds"))

# Save your custom plots:
ggsave(paste0(project_dir,"custom_plot.png"), plot = holder)
```