ch13 clustering.Rmd

---
title: "ch13 clustering"
author: "Bernie Mulvey"
date: "2023-05-15"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.height = 10,fig.width = 7,include = FALSE)
knitr::opts_chunk$set(fig.width=7,fig.height=10)
#### sets tab autocompletion for directories to begin from the directory containing the .Rproj session, instead of the script's directory if different
knitr::opts_knit$set(root.dir = here::here())

library(ggplot2)
library(data.table)
library(Biostrings)
library(gridExtra)
require(colorout)
ColorOut()
options("styler.addins_style_transformer" = "biocthis::bioc_style()")
library(SpatialExperiment)
library(ggspavis)
library(scater) # addPerCellQC
library(nnSVG)
library(BiocParallel)
library(scran)
library(parallel)
library(scCustomize)
library(fasthplus)
# library(STexampleData)

theme_set(theme_bw()+theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16), axis.text.y = element_text(size = 14), axis.title.y = element_text(size =16), plot.title = element_text(size = 20,hjust=0.5), strip.text = element_text(size=18), legend.text = element_text(size=10), legend.title = element_text(size=11,hjust=0.5)))
```

load previous objs
```{r}
# filtered data (cloned into the three list objs below, which also included PCA/UMAP info )
# splc2 <- readRDS("splc2_line201.RDS")
# hvgs, svgs, and combinations
vgs <- fread("Top 1310 ea SVGs and HVGs, union, and n union top SVGs.txt")

# umap and pca objs
umaplist <- readRDS("processed/ch12/umaplist.RDS")
# umdirlist <- readRDS("processed/ch12/umapdirlist.RDS")
pcalist <- readRDS("processed/ch12/pcalist.RDS")
```

Spatially-UNAWARE clustering:
"We can perform clustering by applying standard clustering methods developed for single-cell RNA sequencing data, using molecular features (gene expression). Here, we apply graph-based clustering using the Walktrap method implemented in scran (Lun, McCarthy, and Marioni 2016), applied to the top 50 PCs calculated on the set of top HVGs."

Since we have H, S, and H+S vgs, we can try it w all four. the bigger question will be what to set k.

```{r}
bpparam <- MulticoreParam(workers=8)
glist <- bplapply(pcalist,BPPARAM = bpparam,FUN=function(x){
  scran::buildSNNGraph(x,k=16,use.dimred=reducedDimNames(x))
})

g_walk <- bplapply(glist,igraph::cluster_walktrap,BPPARAM = bpparam)
walkclusts <- lapply(g_walk,FUN=function(x){x$membership})

# output is a vector of length(colData((spe)) with the cluster assignment 
walkclusts[[1]]

## check cluster sizes
lapply(walkclusts,table)
## notice that despite k set to 16 above, the different gene sets give 8-14 clusters. not sure what to make of that. 

# append cluster labels to spes for visualization
names(walkclusts) <- names(pcalist)
i<-1
for (i in c(1:4)){
  colLabels(pcalist[[i]]) <- factor(walkclusts[[i]])
}
rm(i,g_walk,glist)
gc(full=T)
# colLabels gives this column the name label and no way to do otherwise
names(colData(pcalist[[1]]))
# what does it do when you feed a second label in?
colLabels(pcalist[[1]]) <- rep("z",ncol(pcalist[[1]]))
colData(pcalist[[1]])[1:3,32]
## oh, just overwrites it. so hang onto walkclusts so we can rotate through them and other labels as we go. and switch the pcal [[1]] back to the walkclust assignments
colLabels(pcalist[[1]]) <- factor(walkclusts[[1]])
# saveRDS(walkclusts,"data/pcalist_walktrap_clusts.RDS")
```

plot cluster assignments onto spatial data
### note: scCustomize has functions to pull from palette sets of up to 50 colors
### https://samuel-marsh.github.io/scCustomize/articles/Color_Palettes.html
### and see here for storing those palette assignments https://samuel-marsh.github.io/scCustomize/articles/Helpers_and_Utilities.html
```{r}

# ick, the first two colors are black and gray so skip those ones by pulling 2 extra colors and only using colors ≥3
plts <- list()
i<-1
for (i in c(1:length(pcalist))){
  x<-pcalist[[i]]
  plts[[i]] <- ggspavis::plotSpots(pcalist[[i]], annotate = "label", palette = DiscretePalette_scCustomize(num_colors = 2 + length(unique(colData(x)$label)), palette = "polychrome")[3:(length(unique(colData(x)$label)) + 2)])+
    ggtitle(paste0(names(pcalist)[i]))
  rm(x)
}
rm(i)

# pdf("plots/ch13_clusters_walktrap_k16_by_hsvgSets.pdf",height=18,width=4)
do.call("grid.arrange",c(plts,ncol=1))
# dev.off()
rm(plts)

```


### now i'm not really sure what this walktrap business is...i think we need to play with other spatially-aware and unaware clustering algos. let's look through the LC and dlPFC git repos for what some more refined clustering techniques might look like.



## 05/18/23 ##
```{r}
# restore stuff from above + umaps derived from pcalist
pcalist <- readRDS("processed/ch12/pcalist.RDS")
umaplist <- readRDS("processed/ch12/umaplist.RDS")
walkclusts <- readRDS("data/pcalist_walktrap_clusts.RDS")
i<-1
for (i in c(1:length(pcalist))){
  colLabels(pcalist[[i]]) <- factor(walkclusts[[i]])
  colLabels(umaplist[[i]]) <- factor(walkclusts[[i]])
  reducedDimNames(umaplist[[i]])[1] <- "PCA"
  reducedDimNames(umaplist[[i]])[2] <- "UMAP"}

plts <- list()
i<-1
js <- seq(1,by=2,to=2*length(umaplist))
for (i in c(1:length(umaplist))){
  j <- js[i]
  x <- umaplist[[i]]
  plts[[j]] <-
    plotDimRed(
    x,
    type = "PCA",
    annotate = "label",
    palette = DiscretePalette_scCustomize(num_colors = 2 + length(unique(colData(x)$label)), palette = "polychrome")[3:(length(unique(colData(x)$label)) + 2)]
) + aes(shape=factor(colData(x)$donor_id))
    ggtitle(paste0(gsub(
        names(umaplist)[i],
        pattern = "_UMAP", replacement = ""
    )))
  
# in umap space
  plts[[j + 1]] <-
      plotDimRed(
          x,
          type = "UMAP",
          annotate = "label",
          palette = DiscretePalette_scCustomize(num_colors = 2 + length(unique(colData(
              x
          )$label)), palette = "polychrome")[3:(length(unique(colData(x)$label)) + 2)]
      ) +
      ggtitle(names(umaplist)[i])
  rm(x,j)
}

rm(i,js)

# pdf("plots/ch13_clustersInDimRed_walktrap_k16_by_hsvgs.pdf",height=18,width=10)
do.call("grid.arrange",c(plts,nrow=4))
# dev.off()
rm(plts)


#### for reference: though ugly, to add shapes coding explanatory variables:
# plotDimRed(
#     x,
#     type = "PCA",
#     annotate = "label",
#     palette = DiscretePalette_scCustomize(num_colors = 2 + length(unique(colData(x)$label)), palette = "polychrome")[3:(length(unique(colData(x)$label)) + 2)],
#     size = 1
# ) + aes(shape=factor(colData(x)$donor_id)) +
# scale_shape_manual(values=seq(7,to=length(colData(x)$donor_id)))

## the scale shape line grabs as many shapes as necess for the explanatory var, so long as its under 25 levels anyhow
```


### Extending from ch13, part 1: optimal k determination w/ fasthplus
##### adapting from louise H's script for dlPFC spatial: https://github.com/LieberInstitute/spatialDLPFC/blob/main/code/analysis/06_fasthplus/01_fasthplus.R

the 1310 svgs (umap plot 4) look cleanest so we'll use svgs as our feature set here.
```{r}
pcalist <- pcalist[[4]]
gc(full=T)
### we can do 5 at a time locally with a dataset this size, so:
kspan <- c(6:30)
fhrestab <- as.data.frame(matrix(nrow=length(kspan),ncol=2))
colnames(fhrestab) <- c("k","fasthplus")

i<-1
for (i in c(1:length(kspan)/5)){
  ks <- kspan[c((5*i+1):(5*i+5))]
  multiklist <- list(pcalist,pcalist,pcalist,pcalist,pcalist)
  multiklist <- lapply(multiklist,FUN=function(x){
    reducedDimNames(x) <- "PCA"
    x})
  # hpb estimate. t = pre-bootstrap sample size, D = reduced dimensions matrix, L = cluster labels, r = number of bootstrap iterations
  
  # helper functions
  find_t <- function(L, proportion = 0.05) {
        initial_t <- floor(length(L) * proportion)
        smallest_cluster_size <- min(table(L))
        n_labels <- length(unique(L))
        ifelse(smallest_cluster_size > (initial_t / n_labels), initial_t, smallest_cluster_size * n_labels)
  }
  
  bpparam <- MulticoreParam(workers=8)
  register(bpparam)
  
  mk.glist <- bpmapply(function(X,Y){scran::buildSNNGraph(X,k=Y,use.dimred=reducedDimNames(X))}, BPPARAM = bpparam, X=multiklist, Y=ks)
  
  mk.g_walk <- bplapply(mk.glist,igraph::cluster_walktrap,BPPARAM = bpparam)
  rm(mk.glist)
  mk.walkclusts <- lapply(mk.g_walk,FUN=function(x){x$membership})
  rm(mk.g_walk)
  
  # assign as character, not factor -- otherwise we can't drop rows properly in fasth loop below
  multiklist <- bpmapply(X=mk.walkclusts,Y=multiklist,FUN=function(X,Y){
    colLabels(Y) <- factor(X)
    Y})
  # check if all this worked
  lapply(multiklist,FUN=function(x){table(colData(x)$label)})
  # ok
  
  multiklist <- mapply(X=multiklist,Y=ks,FUN=function(X,Y){
    colnames(colData(X))[which(colnames(colData(X))=="label")] <- paste0("label_",Y)
    X
  })
  
  fhres <- mapply(X=multiklist,Y=ks,function(X,Y){
    initial_t <- find_t(L = colData(X)[[paste0("label_",Y)]], proportion = 0.01)
    cluster_prop <- table(colData(X)[[paste0("label_",Y)]]) / ncol(X)
    bad_clusters <- which(cluster_prop < 0.01)
      if (length(bad_clusters) > 0) {
          message("For k: ", Y, " we are dropping small clusters: ", paste(names(bad_clusters), collapse = ", "))
          X <- X[, !(levels(colData(X)[[paste0("label_",Y)]]) %in% as.character(names(bad_clusters)))]
          updated_t <- find_t(colData(X)[[paste0("label_",Y)]], 0.01)
          message("initial t: ", initial_t, "; updated t: ", updated_t)
      }
      else{
        updated_t <- initial_t
      }
      
    set.seed(42)
    hpb(D = reducedDims(X)$PCA, L = colData(X)[[paste0("label_", Y)]], t = updated_t, r = 100)
  })
  
  if(i==1){fhrestab <- cbind(ks,unlist(fhres))}
  else{fhrestab <- rbind(fhrestab,cbind(ks,unlist(fhres)))}
  rm(curks,fhres,multiklist,mk.walkclusts)
  gc(full=T)
}
write.table(fhrestab,"analysis/ch13-clustering/fasthplus_results_1310svgs.txt",sep='\t',quote=F,row.names=F,col.names = T)


```