---
title: "Suryawanshi_2018"
output: html_document
---
```{r}
# Load the libraries
library(SingleCellExperiment)
library(Seurat)
remotes::install_github("mojaveazure/seurat-disk")
library(SeuratDisk)
library(tidyverse)
library(RColorBrewer)
library(scDblFinder)
library(cowplot)
```
The very first thing to do is load the data.
```{r}
# Use the connect() function to read in the loom file
loom_obj <- Connect(filename = "first-trimester-decidua-human-placenta-10XV2.loom", mode = 'r')

seurat_loom <- as.Seurat(loom_obj)
saveRDS(seurat_loom, "seuratObj.rds")
```

Let's perform basic QC steps
```{r}
View(seurat_loom@meta.data)
median(seurat_loom@meta.data$nCount_RNA)
median(seurat_loom@meta.data$nFeature_RNA)

seurat_loom$percent.mt <- PercentageFeatureSet(seurat_loom, pattern = "^MT-")

# nCount_RNA is the number of UMI counts in a cell
hist(seurat_loom$nCount_RNA)

# nFeature_RNA is the number of different genes that had any reads
hist(seurat_loom$nFeature_RNA)

# percent.mt is the percent mitochondrial reads
hist(seurat_loom$percent.mt)

# Make a violin plot of the QC columns
plt <- VlnPlot(seurat_loom, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
               ncol = 3) 
plt
ggsave(filename = "QC.png", plot = plt, width = 7, height = 3.5)

# Make scatter plots to compare
plot1 <- FeatureScatter(seurat_loom, 
                        feature1 = "nCount_RNA", 
                        feature2 = "percent.mt")
plot1
plot2 <- FeatureScatter(seurat_loom, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot2

plot_grid(plot1, plot2, nrow = 2)


# Subset the seurat_loom object by reasonable thresholds
dim(seurat_loom)
seurat_loom <- subset(seurat_loom, subset = nFeature_RNA > 200 &
                 nFeature_RNA < 6500 & percent.mt < 15)

dim(seurat_loom)

saveRDS(seurat_loom, "Filtered_seuratObj.rds")
```

Now we proceed to normalization of the data
```{r}
# Log-transform the counts
seurat_loom <- NormalizeData(seurat_loom)

# Find Variable Features
seurat_loom <- FindVariableFeatures(seurat_loom)

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

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(seurat_loom)
plot1 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1

# Scale the data
seurat_loom <- ScaleData(seurat_loom)

# Run PCA
seurat_loom <- RunPCA(seurat_loom)

# Plot the PCA
DimPlot(seurat_loom, reduction = "pca")

# Plot the loadings
VizDimLoadings(seurat_loom, dims = 1:2, reduction = "pca")

# Choose the number of principle components to keep
ElbowPlot(seurat_loom)

```