97 lines (73 with data), 2.7 kB

---
title: "R Notebook"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. 

```{r}
# Create copy of spe for exploration
spe_sub_copy <- spe_sub
```

### Exploring Zero-Cell and Single-Cell Spots
```{r}
# distribution of cells per spot
tbl_cells_per_spot[1:13]

# as proportions
prop_cells_per_spot <- round(tbl_cells_per_spot / sum(tbl_cells_per_spot), 2)
prop_cells_per_spot[1:13]
```

## Feature Selection: Identify Highly Variable Genes (HVGs)

```{r}
# remove mitochondrial genes
spe_sub_copy <- spe_sub_copy[!is_mito,]
dim(spe_sub_copy)

library(scran)

# fit mean-variance relationship
dec <- modelGeneVar(spe_sub_copy)

# visualize mean-variance relationship
fit <- metadata(dec)
plot(fit$mean, fit$var, 
     xlab = "mean of log-expression", ylab = "variance of log-expression")
curve(fit$trend(x), col = "dodgerblue", add = TRUE, lwd = 2)

# select top HVGs
top_hvgs <- getTopHVGs(dec, prop = 0.1)
length(top_hvgs)
```

```{r}
# compute PCA
set.seed(123)
spe_sub_copy <- runPCA(spe_sub_copy, subset_row = top_hvgs)
# compute UMAP on top 50 PCs
set.seed(123)
spe_sub_copy <- runUMAP(spe_sub_copy, dimred = "PCA")
# update column names
colnames(reducedDim(spe_sub_copy, "UMAP")) <- paste0("UMAP", 1:2)

# plot top 2 PCA dimensions
plotDimRed(spe_sub_copy, type = "PCA")

# plot top 2 UMAP dimensions
plotDimRed(spe_sub_copy, type = "UMAP")
```

```{r}
# graph-based clustering
set.seed(123)
k <- 12
g <- buildSNNGraph(spe_sub_copy, k = k, use.dimred = "PCA")
g_walk <- igraph::cluster_walktrap(g)
```

```{r}
clus <- igraph::cut_at(g_walk, no = 6)
table(clus)
# store cluster labels in column 'label' in colData
colLabels(spe) <- factor(clus)
```

```{r}
# plot clusters in PCA reduced dimensions
plotDimRed(spe_sub_copy,type = "PCA", 
           annotate = "label", palette = "libd_layer_colors")

# plot clusters in UMAP reduced dimensions
plotDimRed(spe_sub_copy, type = "UMAP", 
           annotate = "label", palette = "libd_layer_colors")
```


Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Cmd+Shift+K* to preview the HTML file). 

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.