a b/HVGs_DimReduce.Rmd 

  1
  
---





  
  

  2
  
title: "R Notebook"





  
  

  3
  
output: html_notebook





  
  

  4
  
---





  
  

  5
  





  
  

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





  
  

  7
  





  
  

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





  
  

  9
  





  
  

  10
  
```{r}





  
  

  11
  
# Create copy of spe for exploration





  
  

  12
  
spe_sub_copy <- spe_sub





  
  

  13
  
```





  
  

  14
  





  
  

  15
  
### Exploring Zero-Cell and Single-Cell Spots





  
  

  16
  
```{r}





  
  

  17
  
# distribution of cells per spot





  
  

  18
  
tbl_cells_per_spot[1:13]





  
  

  19
  





  
  

  20
  
# as proportions





  
  

  21
  
prop_cells_per_spot <- round(tbl_cells_per_spot / sum(tbl_cells_per_spot), 2)





  
  

  22
  
prop_cells_per_spot[1:13]





  
  

  23
  
```





  
  

  24
  





  
  

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





  
  

  26
  





  
  

  27
  
```{r}





  
  

  28
  
# remove mitochondrial genes





  
  

  29
  
spe_sub_copy <- spe_sub_copy[!is_mito,]





  
  

  30
  
dim(spe_sub_copy)





  
  

  31
  





  
  

  32
  
library(scran)





  
  

  33
  





  
  

  34
  
# fit mean-variance relationship





  
  

  35
  
dec <- modelGeneVar(spe_sub_copy)





  
  

  36
  





  
  

  37
  
# visualize mean-variance relationship





  
  

  38
  
fit <- metadata(dec)





  
  

  39
  
plot(fit$mean, fit$var, 





  
  

  40
  
     xlab = "mean of log-expression", ylab = "variance of log-expression")





  
  

  41
  
curve(fit$trend(x), col = "dodgerblue", add = TRUE, lwd = 2)





  
  

  42
  





  
  

  43
  
# select top HVGs





  
  

  44
  
top_hvgs <- getTopHVGs(dec, prop = 0.1)





  
  

  45
  
length(top_hvgs)





  
  

  46
  
```





  
  

  47
  





  
  

  48
  
```{r}





  
  

  49
  
# compute PCA





  
  

  50
  
set.seed(123)





  
  

  51
  
spe_sub_copy <- runPCA(spe_sub_copy, subset_row = top_hvgs)





  
  

  52
  
# compute UMAP on top 50 PCs





  
  

  53
  
set.seed(123)





  
  

  54
  
spe_sub_copy <- runUMAP(spe_sub_copy, dimred = "PCA")





  
  

  55
  
# update column names





  
  

  56
  
colnames(reducedDim(spe_sub_copy, "UMAP")) <- paste0("UMAP", 1:2)





  
  

  57
  





  
  

  58
  
# plot top 2 PCA dimensions





  
  

  59
  
plotDimRed(spe_sub_copy, type = "PCA")





  
  

  60
  





  
  

  61
  
# plot top 2 UMAP dimensions





  
  

  62
  
plotDimRed(spe_sub_copy, type = "UMAP")





  
  

  63
  
```





  
  

  64
  





  
  

  65
  
```{r}





  
  

  66
  
# graph-based clustering





  
  

  67
  
set.seed(123)





  
  

  68
  
k <- 12





  
  

  69
  
g <- buildSNNGraph(spe_sub_copy, k = k, use.dimred = "PCA")





  
  

  70
  
g_walk <- igraph::cluster_walktrap(g)





  
  

  71
  
```





  
  

  72
  





  
  

  73
  
```{r}





  
  

  74
  
clus <- igraph::cut_at(g_walk, no = 6)





  
  

  75
  
table(clus)





  
  

  76
  
# store cluster labels in column 'label' in colData





  
  

  77
  
colLabels(spe) <- factor(clus)





  
  

  78
  
```





  
  

  79
  





  
  

  80
  
```{r}





  
  

  81
  
# plot clusters in PCA reduced dimensions





  
  

  82
  
plotDimRed(spe_sub_copy,type = "PCA", 





  
  

  83
  
           annotate = "label", palette = "libd_layer_colors")





  
  

  84
  





  
  

  85
  
# plot clusters in UMAP reduced dimensions





  
  

  86
  
plotDimRed(spe_sub_copy, type = "UMAP", 





  
  

  87
  
           annotate = "label", palette = "libd_layer_colors")





  
  

  88
  
```





  
  

  89
  





  
  

  90
  





  
  

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





  
  

  92
  





  
  

  93
  
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). 





  
  

  94
  





  
  

  95
  
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.





  
  

  96