Additionally, we can use regression to remove any unwanted sources of variation from the dataset, such as `cell cycle`, `sequencing depth`, `percent mitocondria`. This is achieved by doing a generalized linear regression using these parameters as covariates in the model. Then the residuals of the model are taken as the "regressed data". Although perhaps not in the best way, batch effect regression can also be done here. ```{r} alldata <- ScaleData(alldata, vars.to.regress = c("percent_mito", "nFeature_RNA"), assay = "RNA") ``` ## PCA *** Performing PCA has many useful applications and interpretations, which much depends on the data used. In the case of life sciences, we want to segregate samples based on gene expression patterns in the data. ```{r} alldata <- RunPCA(alldata, npcs = 50, verbose = F) ``` We then plot the first principal components. ```{r, fig.asp=.28} plot_grid(ncol = 3, DimPlot(alldata, reduction = "pca", group.by = "orig.ident",dims = 1:2), DimPlot(alldata, reduction = "pca", group.by = "orig.ident",dims = 3:4), DimPlot(alldata, reduction = "pca", group.by = "orig.ident",dims = 5:6) ) ``` To identify which genes (Seurat) or metadata paramters (Scater/Scran) contribute the most to each PC, one can retreive the loading matrix information. Unfortunatelly this is not implemented in Scater/Scran, so you will need to compute PCA using `logcounts`. ```{r,fig.asp=.5} VizDimLoadings(alldata, dims = 1:5, reduction = "pca",ncol = 5,balanced = T) ``` We can also plot the amount of variance explained by each PC. ```{r,fig.asp=.3} ElbowPlot(alldata, reduction = "pca",ndims = 50) ``` Based on this plot, we can see that the top 8 PCs retain a lot of information, while other PCs contain pregressivelly less. However, it is still advisable to use more PCs since they might contain informaktion about rare cell types (such as platelets and DCs in this dataset) ## tSNE *** We can now run [BH-tSNE](https://arxiv.org/abs/1301.3342). ```{r,fig.asp=1} alldata <- RunTSNE(alldata, reduction = "pca", dims = 1:30, perplexity=30, max_iter=1000, theta=0.5, eta=200, num_threads=0 ) #see ?Rtsne and ?RunTSNE for more info ``` We can now plot the tSNE colored per dataset. We can clearly see the effect of batches present in the dataset. ```{r,fig.asp=.28} plot_grid(ncol = 3,DimPlot(alldata, reduction = "tsne", group.by = "orig.ident")) ``` *** ## UMAP *** We can now run [UMAP](https://arxiv.org/abs/1802.03426) for cell embeddings. ```{r} alldata <- RunUMAP(alldata, reduction = "pca", dims = 1:30, n.components=2, n.neighbors=30, n.epochs=200, min.dist=0.3, learning.rate=1, spread=1 ) #see ?RunUMAP for more info ```

**Your turn**
We have now done Variable gene selection, PCA and UMAP with the settings we chose. Test a few different ways of selecting variable genes, number of PCs for UMAP and check how it influences your embedding.

### TASK: Test settings
Another usefullness of UMAP is that it is not limitted by the number of dimensions the data cen be reduced into (unlike tSNE). We can simply reduce the dimentions altering the `n.components` parameter.
```{r}
# we can add in additional reductions, by defulat they are named "pca", "umap", "tsne" etc. But we can specify alternative names with reduction.name
alldata <- RunUMAP(alldata, reduction.name = "UMAP10_on_PCA",
reduction = "pca",
dims = 1:30,
n.components=10,
n.neighbors=30,
n.epochs=200,
min.dist=0.3,
learning.rate=1,
spread=1 )
#see ?RunUMAP for more info
```
We can now plot the UMAP colored per dataset. Although less distinct as in the tSNE, we still see quite an effect of the different batches in the data.
```{r,fig.asp=.28}
plot_grid(ncol = 3,
DimPlot(alldata, reduction = "umap", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_PCA"),
DimPlot(alldata, reduction = "UMAP10_on_PCA", group.by = "orig.ident",dims = 1:2)+ ggplot2::ggtitle(label ="UMAP10_on_PCA"),
DimPlot(alldata, reduction = "UMAP10_on_PCA", group.by = "orig.ident",dims = 3:4)+ ggplot2::ggtitle(label ="UMAP10_on_PCA")
)
```
We can now plot PCA, UMAP and tSNE side by side for comparison. Here, we can conclude that our dataset contains a batch effect that needs to be corrected before proceeding to clustering and differential gene expression analysis.
```{r,fig.asp=.28}
plot_grid(ncol = 3,
DimPlot(alldata, reduction = "pca", group.by = "orig.ident"),
DimPlot(alldata, reduction = "tsne", group.by = "orig.ident"),
DimPlot(alldata, reduction = "umap", group.by = "orig.ident")
)
```
## Using ScaledData and graphs for DR
***
Althought running a sencond dimmensionality reduction (i.e tSNE or UMAP) on PCA would be a standard approach (because it allows higher computation efficiency), the options are actually limiteless. Below we will show a couple of other common options such as running directly on the scaled data (which was used for PCA) or on a graph built from scaled data. We will show from now on only UMAP, but the same applies for tSNE.
### Using ScaledData for UMAP
To run tSNE or UMAP on the scaled data, one firts needs to select the number of variables to use. This is because including dimentions that do contribute to the separation of your cell types will in the end mask those differences. Another reason for it is because running with all genes/features also will take longer or might be computationally unfeasible. Therefore we will use the scaled data of the highly variable genes.
```{r}
alldata <- RunUMAP(alldata, reduction.name = "UMAP_on_ScaleData",
features = alldata@assays$RNA@var.features,
assay = "RNA",
n.components=2,
n.neighbors=30,
n.epochs=200,
min.dist=0.3,
learning.rate=1,
spread=1 )
```
### Using a Graph for UMAP
To run tSNE or UMAP on the a graph, we first need to build a graph from the data. In fact, both tSNE and UMAP first build a graph from the data using a specified distance metrix and then optimize the embedding. Since a graph is just a matrix containing distances from cell to cell and as such, you can run either UMAP or tSNE using any other distance metric desired. Euclidean and Correlation are ususally the most commonly used.
```{r}
#Build Graph
alldata <- FindNeighbors(alldata,
reduction = "pca",
graph.name = "SNN",
assay = "RNA",
k.param = 20,
features = alldata@assays$RNA@var.features)
#Run UMAP on a graph
alldata <- RunUMAP(alldata, reduction.name = "UMAP_on_Graph",
graph = "SNN",
assay = "RNA" )
```
We can now plot the UMAP comparing both on PCA vs ScaledSata vs Graph.
```{r, fig.asp=.28}
p1 <- DimPlot(alldata, reduction = "umap", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_PCA")
p2 <- DimPlot(alldata, reduction = "UMAP_on_ScaleData", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_ScaleData")
p3 <- DimPlot(alldata, reduction = "UMAP_on_Graph", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_Graph")
leg <- get_legend(p1)
gridExtra::grid.arrange(
gridExtra::arrangeGrob(
p1 + NoLegend() + NoAxes(),
p2 + NoLegend() + NoAxes(),
p3 + NoLegend() + NoAxes(),
leg,nrow=2),
ncol=1,widths=c(1)
)
```
## Ploting genes of interest
***
Let's plot some marker genes for different celltypes onto the embedding. Some genes are:
Markers | Cell Type
--- | ---
CD3E | T cells
CD3E CD4 | CD4+ T cells
CD3E CD8A | CD8+ T cells
GNLY, NKG7 | NK cells
MS4A1 | B cells
CD14, LYZ, CST3, MS4A7 | CD14+ Monocytes
FCGR3A, LYZ, CST3, MS4A7 | FCGR3A+ Monocytes
FCER1A, CST3 | DCs
```{r,fig.asp=.9}
myfeatures <- c("CD3E","CD4","CD8A","NKG7","GNLY","MS4A1","CD14","LYZ","MS4A7","FCGR3A","CST3","FCER1A")
FeaturePlot(alldata, reduction = "umap",dims = 1:2,
features = myfeatures,ncol = 3,order = T)
for(i in myfeatures){
}
p1 <- DimPlot(alldata, reduction = "umap", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_PCA")
p2 <- DimPlot(alldata, reduction = "UMAP_on_ScaleData", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_ScaleData")
p3 <- DimPlot(alldata, reduction = "UMAP_on_Graph", group.by = "orig.ident")+ ggplot2::ggtitle(label ="UMAP_on_Graph")
leg <- get_legend(p1)
gridExtra::grid.arrange(
gridExtra::arrangeGrob(
p1 + NoLegend() + NoAxes(),
p2 + NoLegend() + NoAxes(),
p3 + NoLegend() + NoAxes(),
leg,nrow=2),
ncol=1,widths=c(1)
)
```
We can finally save the object for use in future steps.
```{r}
saveRDS(alldata,"data/results/covid_qc_dr.rds")
```
**Your turn**
Select some of your dimensionality reductions and plot some of the QC stats that were calculated in the previous lab. Can you see if some of the separation in your data is driven by quality of the cells?

### Session Info
***
```{r}
sessionInfo()
```