3D visualization of STitch3D’s generated virtual slice

rm(list = ls())
options(warn=-1)
library(rgl)
library(misc3d)
knitr::knit_hooks$set(webgl = hook_webgl)

Load plot3D functions.

source("./plot3D_func.R")

Load cell-type proportions and 3D coordinates of spots

directory = "./res_mouse_brain"
file_list <- list.files(path = directory)
n_slice <- sum(unlist(lapply(file_list, function(x){startsWith(x, "prop_slice")})))
#cell-type proportions
for (i in (0:(n_slice-1))){
    prop <- read.table(paste0(directory, "/prop_slice", i, ".csv"), sep=",", header=TRUE)
    if (i == 0){
        prop_all <- prop
    }else{
        prop_all <- rbind(prop_all, prop)
    }
}
colnames(prop_all)[1] <- "spot"
#3D coordinates
coor_3d <- read.table(paste0(directory, "/3D_coordinates.csv"), sep=",", header=TRUE)
colnames(coor_3d)[1] <- "spot"

spots.table <- merge(coor_3d, prop_all, by=c("spot"))

alpha_threshold <- 0.2
alpha_background <- 0.02

Generate the virtual slice

a <- 1
b <- 0
c <- 0
d <- -28
n <- c(a,b,c)/sqrt(a^2+b^2+c^3)
scale <- 1.1

dist <- abs(a*spots.table$x + b*spots.table$y + c*spots.table$z + d) / sqrt(a^2+b^2+c^2)
p <- - (a*spots.table$x + b*spots.table$y + c*spots.table$z + d) / (a^2+b^2+c^2)
spots.table$dist <- dist
spots.table$p <- p

STitch3D generates the virtual slice.

Distributions of four hippocampal neuron types

celltypes <- c("Ext_Hpc_CA1", "Ext_Hpc_CA2", "Ext_Hpc_CA3", "Ext_Hpc_DG1")
celltype_colors <- c("#d62728", "#FFFF00", "#1f77b4", "#2ca02c")
um <- c(0.98223096, -0.05770408, -0.1785820, 0,
        0.02726252, -0.89759272, 0.4399813, 0,
        -0.18568322, -0.43703181, -0.8800704, 0,
        0, 0, 0, 1) #set the initial view of the 3D plot

prop <- spots.table[, celltypes]
prop$max_prop <- apply(prop, 1, max)
prop$max_celltype <- apply(prop, 1, function(x){celltypes[which.max(x)]})
prop$color <- "gray"
prop$alpha <- prop$max_prop

#for others
prop$alpha[prop$max_prop <= alpha_threshold] <- alpha_background
#for target cell types
for (i in 1:length(celltypes)){
      prop$color[(prop$max_prop > alpha_threshold) & (as.vector(prop$max_celltype) == celltypes[i])] <- celltype_colors[i]
}

#3D plot
open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   1

par3d(persp)

## NULL

view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))
spheres3d(spots.table$x, spots.table$y, spots.table$z,
          col = prop$color, radius=0.5, 
          alpha = prop$alpha)
decorate3d()

rgl.planes(a=a, b=b, c=c, d=d, alpha=1, color="gray", lit=FALSE, smooth=FALSE, depth_test='always',
           xlim = c(min(spots.table$x), max(spots.table$x)), 
           ylim = c(min(spots.table$y), max(spots.table$y)), 
           zlim = c(min(spots.table$z), max(spots.table$z)))

um <- c(0, 0, -1, 0,
        0, -1, 0, 0,
        -1, 0, 0, 0,
        0, 0, 0, 1)
subdf <- spots.table[spots.table$dist <= scale, ]
subprop <- prop[spots.table$dist <= scale, ]
subprop$alpha[subprop$color == "gray"] <- 0.1

open3d(windowRect = c(0, 0, 720, 1080))

## glX 
##   3

par3d(persp)

## NULL

view3d(userMatrix=matrix(um, byrow=TRUE, nrow=4))

spheres3d(subdf$x, subdf$y, subdf$z,
          col = subprop$color, radius=0.5, 
          alpha = subprop$alpha)

Distributions of two thalamic neuron types

celltypes <- c("Ext_Thal_1", "Ext_Thal_2")
celltype_colors <- c("orange", "#008B8B")
um <- c(0.98223096, -0.05770408, -0.1785820, 0,
        0.02726252, -0.89759272, 0.4399813, 0,
        -0.18568322, -0.43703181, -0.8800704, 0,
        0, 0, 0, 1) #set the initial view of the 3D plot

prop <- spots.table[, celltypes]
prop$max_prop <- apply(prop, 1, max)
prop$max_celltype <- apply(prop, 1, function(x){celltypes[which.max(x)]})
prop$color <- "gray"
prop$alpha <- prop$max_prop

#for others
prop$alpha[prop$max_prop <= alpha_threshold] <- alpha_background
#for target cell types
for (i in 1:length(celltypes)){
      prop$color[(prop$max_prop > alpha_threshold) & (as.vector(prop$max_celltype) == celltypes[i])] <- celltype_colors[i]
}

#3D plot
open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   4

par3d(persp)

## NULL

view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))
spheres3d(spots.table$x, spots.table$y, spots.table$z,
          col = prop$color, radius=0.5, 
          alpha = prop$alpha)
decorate3d()

rgl.planes(a=a, b=b, c=c, d=d, alpha=1, color="gray", lit=FALSE, smooth=FALSE, depth_test='always',
           xlim = c(min(spots.table$x), max(spots.table$x)), 
           ylim = c(min(spots.table$y), max(spots.table$y)), 
           zlim = c(min(spots.table$z), max(spots.table$z)))

um <- c(0, 0, -1, 0,
        0, -1, 0, 0,
        -1, 0, 0, 0,
        0, 0, 0, 1)
subdf <- spots.table[spots.table$dist <= scale, ]
subprop <- prop[spots.table$dist <= scale, ]
subprop$alpha[subprop$color == "gray"] <- 0.1

open3d(windowRect = c(0, 0, 720, 1080))

## glX 
##   6

par3d(persp)

## NULL

view3d(userMatrix=matrix(um, byrow=TRUE, nrow=4))

spheres3d(subdf$x, subdf$y, subdf$z,
          col = subprop$color, radius=0.5, 
          alpha = subprop$alpha)