Graph Based Clustering

##################
# LOAD LIBRARIES #
##################
#setwd("/Users/vilkal/work/Brolidens_work/Projects/Gabriella_repo/reports/rmarkdown/manuscript")
suppressWarnings({suppressMessages({suppressPackageStartupMessages({
  library(tidyverse)
  library(igraph)
  library(rafalib)
  #remotes::install_github("czarnewski/niceRplots",upgrade = "always",force = T)
  library(niceRplots)
})  })  })

#########
# PATHS #
#########
input_dir <- "../results/03_normalize_data_output/"
result_dir <- "../results/04_clustering_output/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }

#############
# LODA DATA #
#############
datasets_all_samples <- readRDS(paste0(input_dir, "datasets_all_samples.RDS"))
# datasets_all_samples <- readRDS("../../Gabriella_repo/results/datasets_all_samples.RDS") %>%
#   set_names(., c("Tissue_RNAseq_V3_normalized", "ASV_Tissue_normalized", "ASV_Luminal_normalized", "ASV_CVL_V2_normalized_NOT_batch_corrected"))
metadata <- read.csv(paste0("../data/", "metadata.csv"),row.names = 1)


#########################
# SELECT SAMPLES TO USE #
#########################
# datasets_all_samples <- datasets_all_samples %>%
#   map(., ~as_tibble(.x, rownames = NA )) %>%
#   map(., ~select(.x, any_of(metadata$ID)))


#################
# COLOR PALETTS #
#################
pal <- c( "#0072B2", "#009E73","#D55E00", "#CC79A7", "#E69F00", "#999999")
taxa_pal <- c(RColorBrewer::brewer.pal(8,"Pastel2"),RColorBrewer::brewer.pal(8,"Pastel1"),"grey90")
bact_pal <- c('#88CCEE', '#44AA99', '#117733', '#332288', '#DDCC77', '#999933','#CC6677', '#882255', '#AA4499', '#DDDDDD')

###################################
# STUDY GROUPS GRAPH CONSTRUCTION #
####################################
all_microbiome <- datasets_all_samples[["ASV_Luminal_normalized"]]
all_microbiome <- all_microbiome[ rowSums(all_microbiome>0) >= 3 ,  ]
dim(all_microbiome)

## [1] 111 108

zscore <- t(apply(all_microbiome,1,function(i){scale(i,T,T)}))
zscore[ zscore > 4 ] <- 4
zscore[ zscore < -4 ] <- -4
rownames(zscore) <- rownames(all_microbiome); colnames(zscore) <- colnames(all_microbiome)

mypar(mar=c(2,15,2,2))
top_var <- apply(all_microbiome,1,var)
top_var_t <- sort(top_var,decreasing = T)[100]
is_top_var <- top_var > top_var_t 
barplot(sort(top_var[is_top_var],decreasing = T)[30:1],horiz = T,las=1)

#Compute sample-wise correlations
cors <- cor(all_microbiome[,], method = "pearson")
adj <- (1-cors)/2
# pheatmap::pheatmap(adj,border_color = NA,color = colorRampPalette(c("grey90","grey70","grey70","grey70","orange3","firebrick"))(90) )

# cors <- cor(all_microbiome, method = "pearson")
# adj <- (1-cors)/2
# pheatmap::pheatmap(adj,border_color = NA,color = colorRampPalette(c("grey90","grey70","grey70","grey70","orange3","firebrick"))(90) )

#Graph construction  
k <- 12 # k for KNN
p <- 0.1 # prunning threshold of SNN 0.14

#Compute KNN graph from correlation matrix
knn <- RANN::nn2(adj,k = k)
nn <- matrix(0, nrow(knn$nn.idx), nrow(knn$nn.idx))
for(idx in 1:nrow(knn$nn.idx) ){  nn[idx,knn$nn.idx[idx,]] <- 1 }

#Computing SNN matrix
N <- matrix(1,ncol = ncol(nn),nrow = nrow(nn))
P <- t(nn) %*% nn
B <- N - nn
C <- t(B) %*% B
Q <- nrow(nn)*N - C
J <- P/Q
J[is.nan(J)] <- 0
J[J < p] <- 0
colnames(J) <- colnames(all_microbiome)
rownames(J) <- colnames(all_microbiome)

#Transforming SNN into an igraph object
g <- graph_from_adjacency_matrix(J,mode = "undirected",diag = F,weighted = T)
write.csv(J,paste0(result_dir,"participant_SNN_graph.csv"),row.names = T)

#####################
# PLOT STUDY GROUPS #
#####################
mypar(3,3,mar=c(0,0,2,0))

set.seed(1)
set.seed(1)
l <- layout_with_fr(g,niter=3000,start.temp=30)

cl <- igraph::cluster_louvain(g)
# cl <- igraph::walktrap.community(g)

metadata$Louvain_clusters <- cl$membership
metadata$joint_clustering <- metadata$Luminal_gr
metadata$joint_clustering <- factor(metadata$joint_clustering) 
set.seed(1)
U <- uwot::umap(adj,n_neighbors = 30)
h <- hclust(as.dist(adj),method = "ward.D2")
# metadata$joint_clustering <- cutree(h,k = 5)

plot( g , vertex.label.cex=0.000001 , vertex.color = pal[factor(cl$membership)] , 
      edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=l,
      vertex.size=10,main="FR embedding\nlouvain",
      edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )

plot( g , vertex.label.cex=0.000001 , vertex.color = pal[factor(cl$membership)] , 
      edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=U,
      vertex.size=10,main="umap",
      edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )
empty_plot()
legend("topleft",title = "Participant groups\nLouvain_gr",
       legend = levels(factor(cl$membership)),
       bty = "n",pch = 21,pt.bg = pal,pt.cex = 1.2,xpd=T)


plot( g , vertex.label.cex=0.000001 , vertex.color = pal[factor(metadata$Luminal_gr)] , 
      edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=l,
      vertex.size=10,main="FR embedding",
      edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )
plot( g , vertex.label.cex=0.000001 , vertex.color = pal[factor(metadata$Luminal_gr)] , 
      edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=U,
      vertex.size=10,main="umap",
      edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )
empty_plot()
legend("topleft",title = "Participant groups\nLuminal_gr",
       legend = levels(factor(metadata$Luminal_gr)),
       bty = "n",pch = 21,pt.bg = pal,pt.cex = 1.2,xpd=T)
# title("Louvain clusters",cex.main=1,font.main=1)


plot( g , vertex.label.cex=0.000001 , vertex.color = pal[factor(metadata$Tissue_gr)] , 
      edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=l,
      vertex.size=10,main="FR embedding",
      edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )
plot( g , vertex.label.cex=0.000001 , vertex.color = pal[factor(metadata$Tissue_gr)] , 
      edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=U,
      vertex.size=10,main="umap",
      edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )
empty_plot()
legend("topleft",title = "Participant groups\nTissue_gr",
       legend = levels(factor(metadata$Tissue_gr)),
       bty = "n",pch = 21,pt.bg = pal,pt.cex = 1.2,xpd=T)

# title("Louvain clusters",cex.main=1,font.main=1)

#######################################
# BACTERIAL GROUPS GRAPH CONSTRUCTION #
#######################################
all_microbiome <- cbind(#datasets_all_samples[["ASV_Tissue_normalized"]],
                        datasets_all_samples[["ASV_Luminal_normalized"]]
                        )

cors <- cor(t(all_microbiome[rowSums(all_microbiome > 0)>=10,]))

#Graph construction 
k_nearest_neighbours <- 5 # k for KNN
snn_threshold <- 0.1 # prunning threshold of SNN

#Compute kNN graph
knn <- RANN::nn2( (1-cors)/2 ,k = k_nearest_neighbours)
nn <- matrix(0, nrow(knn$nn.idx), nrow(knn$nn.idx))
for(idx in 1:nrow(knn$nn.idx) ){  nn[idx,knn$nn.idx[idx,]] <- 1 }

#Calculate SNN and prune based on Jaccard index
N <- matrix(1,ncol = ncol(nn),nrow = nrow(nn))
P <- t(nn) %*% nn
B <- N - nn
C <- t(B) %*% B
Q <- nrow(nn)*N - C
J <- P/Q
J[is.nan(J)] <- 0
J[J < snn_threshold] <- 0
rownames(J) <- rownames(cors)
colnames(J) <- colnames(cors)

#Build a graph object, perform clustering and define the force-directed layout
g <- graph_from_adjacency_matrix(J,mode = "undirected",diag = F,weighted = T)
write.csv(J,paste0(result_dir,"bacteria_SNN_graph.csv"),row.names = T)

cl <- igraph::cluster_louvain(g)
# cl <- igraph::walktrap.community(g)
Bact.com <- tibble(Taxa=rownames(J), Bact_com=cl$membership)
write.csv(Bact.com,paste0(result_dir,"bacterial_communities.csv"),row.names = T)

set.seed(1)
l <- layout_with_fr(g,niter=3000,start.temp=1)

#Define layout with UMAP
mU <- uwot::umap((1-cors)/2, n_neighbors = 10)

##############################
# PLOT BACTERIAL COMMUNITIES #
##############################
mypar(3,3,mar=c(0,0,2,0))
plot( g , vertex.label.cex=0.000001 , vertex.color = bact_pal[factor(cl$membership)] , 
    edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=l,
    vertex.size=10, main="graph layout FR",
    edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )
plot( g , vertex.label.cex=0.000001 , vertex.color = bact_pal[factor(cl$membership)] , 
    edge.width=  ( E(g)$weight / max(E(g)$weight)) ,layout=mU,
    vertex.size=10, main="UMAP",
    edge.color=colorRampPalette(c("grey95","black"))(90) [ round( E(g)$weight / max(E(g)$weight) * 88 )+1 ] )

par(mar=c(0,0,0,0))
empty_plot()
legend("topleft",title = "bact. communities",
     legend = levels(factor(cl$membership)),
     bty = "n",pch = 21,pt.bg = pal,cex = .8,pt.cex = .9,xpd=T)

Figure Bacterial community embedding of 12-SNN graph clustered using Louvain community detection algorithm based on the normalized bacterial abundances across all microbiome datasets (including samples not present in all datasets).