# Set up

To analyze single cell/nucleus data you will most likely need specialized
software (aka. "packages"). In this exercise we will work with one of the most 
popular analysis packages for the R environment called "Seurat". Since we don't 
have time for a full dive into all details of the analysis process, we recommend 
you browse their [vignettes and tutorials](<https://satijalab.org/seurat/articles/get_started_v5_new>). 

For this tutorial, we have preinstalled the packages we need and stored them in the 
`/common/tutorials/bioinformatics/single_cell/software/r-rich4.6_lib/` 
folder in our systems, so you just need to tell R where to find them with the 
`.libPaths` function and make them available for your current session with the 
`library` function: 


```{r libs, warning=FALSE, message=FALSE}

.libPaths() #current software paths 

tutorial_packages <- "/common/tutorials/bioinformatics/single_cell/software/r-rich4.6_lib/" 
.libPaths(c(tutorial_packages, .libPaths())) #adding a new folder
.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, eval=FALSE}

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

# If you want to change it:
project_dir <- "/projects/standard/riss/calix015/bioinfo-training/materials/sc_intro"
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 a directory with the "list.files" function:
list.files(data_dir, recursive = T)

```

# 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:

```{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 using minimum filters, note how we reduced the 
number of genes (rows in the count matrix) from 36601 to 20759:

```{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 
of the contents in a variable instead of its entirety.
We will peek into a few key elements:

```{r so_struct}
# Gene names 
rownames(so) %>% head()

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

# Cell 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, sample integration, clustering, and cluster identification.

Be aware that the thresholds and options chosen in these initial steps are the basis of 
all downstream analysis, so make sure to represent your data with all its 
uniqueness and characteristics. **Think carefully about these steps because they are the base for the results you will get down the line.**


# Pre-analyzed object 

Now we will open a dataset that has been preprocessed and ready to explore and 
generate plots of your favorite genes and cell types.

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

Note that this object has more layers with gene expression information, a set of 
variable genes, dimensional reductions and more information in the metadata 
compared to the one created above:

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

The object metadata 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 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 using the first 2 PCs:

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

Often, you'll see another plot that uses another dimensionality reduction
technique called UMAP, these same cells look like this on a UMAP plot:

```{r umap}
DimPlot(so_done_subset)
```

Note how these two figures are representing exactly the same dataset, but they 
look really different! 

A very **common mistake** is to overinterpret the distance between clusters in a 
UMAP plot and assume they show differences between cells. This is not always
the case! These plots are just a "summary picture" that does not capture 
the whole complexity of the data and that depending on the reduction and options
you use, could look very different. 


## Plot gene expression

Violin and dot plots are often used to summarize gene expression levels of your genes of 
interest across different samples or cell types.

Let's make some with the most variable genes in this data set. By default cells 
are grouped by their identity class, which we have defined above with the `Idents`
function, but you can modify this with the `group.by` parameter.

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

# Dotplot of the top 5 by cell type (Ident):
DotPlot(so_done_subset, top5)

# Violin of the most variable gene by cluster:
VlnPlot(so_done_subset, group.by = "seurat_clusters", top5[1])


```


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}
# See the names of the metadata columns:
colnames(so_done_subset@meta.data)

# Using a metadata column and a gene name:
FeaturePlot(so_done_subset,features = c(top5[1], "nCount_SCT"))

```


## Customizing Plots

Here we will customize the plots made before with the default functions.
You can customize your plots changing or adding function parameters (use 
`?function_name` for details).


```{r custom_umap}
# Get rid of the legend and add labels to cells
DimPlot(so_done_subset,label = TRUE, pt.size = 0.5) + NoLegend()
```

You can also save your plot in a variable and reuse the variable to add custom 
layers:

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

```


# Save your work!

Finally, when you are done with your work remember to save your results and 
plots before your session ends. Whatever is not saved to disk will be lost 
as soon as your session expires.

```{r save, eval=FALSE}
# Save any object with the saveRDS function:
saveRDS(so_done_subset, "subset.rds")

# Save plots with the ggsave function:
ggsave("custom_plot.png", plot = plot_holder)
```