## ----------- 3D Plotly -----------
plot.3d.expr <- function(gene, spots.table, st.data, normalize.by.nuclei = FALSE){
spots.table$expr.level <- st.data[intersect(rownames(st.data),rownames(spots.table)),gene]
spots.plot <- spots.table[intersect(rownames(st.data),rownames(spots.table)),]
# #Discarding outside of brain spots
# print(sprintf('%d spots discarded because they are outside the brain contour.',sum(is.na(spots.plot$acronym))))
# spots.plot <- spots.table[!is.na(spots.plot$acronym),]
#
# #Discarding spots for which gene was not sequenced
# print(sprintf('%d spots discarded because the selected gene was not sequenced.',sum(is.na(spots.plot$expr.level))))
# spots.plot <- spots.plot[!is.na(spots.plot$expr.level),]
#
str.title <- ''
if (normalize.by.nuclei == TRUE){
spots.plot$nuclei[spots.plot$nuclei == 0] <- 1
spots.plot$expr.level <- spots.plot$expr.level / spots.plot$nuclei
str.title <- '(Normalized by number of nuclei)'
}
# change.opacity <- 0.1*max(spots.plot$expr.level)
change.opacity <- quantile(spots.plot$expr.level, 0.8)
spots.plot.1 <- spots.plot[spots.plot$expr.level >= change.opacity,]
spots.plot.2 <- spots.plot[spots.plot$expr.level < change.opacity,]
#Plotting
p <- plot_ly(spots.plot.1,x=~ML,y=~AP,z=~DV,type = 'scatter3d',mode = 'markers',name='High expression',
marker=list(size=2, color = ~expr.level, colorscale = c('#FFE1A1', '#683531'), showscale = TRUE,
opacity = 1,cmin = 0,cmax = max(spots.plot$expr.level)),
hoverinfo = 'text',
text = ~sprintf('Spot: %s<br />DV: %0.2f<br />ML: %0.2f<br />AP: %0.2f<br />Region: %s/%s<br />Expression level: %.0f',
rownames(spots.plot.1),DV,ML,AP,full.name.parent,name,expr.level)) %>%
add_markers(x=spots.plot.2$ML,y=spots.plot.2$AP,z=spots.plot.2$DV,name=sprintf('Low expression (<%0.f)',change.opacity),
marker=list(size=2, color = spots.plot.2$expr.level, colorscale = c('#FFE1A1', '#683531'), showscale = TRUE,
opacity = 0.4,cmin = 0,cmax = max(spots.plot$expr.level),name='Low expression') ,
hoverinfo = 'text',
text = sprintf('Spot: %s<br />DV: %0.2f<br />ML: %0.2f<br />AP: %0.2f<br />Region: %s/%s<br />Expression level: %.0f',
rownames(spots.plot.2),spots.plot.2$DV,spots.plot.2$ML,spots.plot.2$AP,spots.plot.2$full.name.parent,spots.plot.2$name,spots.plot.2$expr.level))%>%
layout(title=sprintf("Gene: %s %s",gene,str.title),
scene = list(xaxis = list(title = 'x = Medial-lateral (mm)',range = c(0,5)),
zaxis = list(title = 'z = Dorsal-ventral (mm)',range = c(-7.9,0)),
yaxis = list(title = 'y = Anterior-posterior (mm)',range = c(-5.9,3))))
p
return(p)
}
## ----------- 2D Plotly -----------
#Plotting expression levels in the atlas
plot.2d.atlas.expr <- function(gene, spots.table, st.data, AP = NULL, slice_index = NULL, min.expr = 0, cex = 1.2, new.window = TRUE, text.plot = TRUE, log.scale = FALSE, normalized.scale = FALSE, win.dim = c(14,14), os.windows = FALSE, spots.to.color = NULL, spots.to.color.2 = NULL, cluster.name.1 = NULL, cluster.name.2 = NULL, X11.type = NULL){
#gene(str): gene to plot
#spots.table(df)
#st.data(df)
#AP(float): AP coordinates. If does not match any sections, takes the closest. If two values, plot all sections in-between.
#slice_index(str): If specifcied, replaces AP coordinates
#min.expr(float): Expression threhsold to plot spot. If between 0 and 1, defines a percentile. Otherwise, defines a number of reads (not log normalized)
#cex(float): dimension of text
#new.window(bool): force new window
#text.plot(bool): show the legend and title
#log.scale(bool): use a logarithmic scale instead
#normalized.scale(bool): all expression are a value between 0 and 100 (maximum expression across the whole brain)
#win.dim(numeric vector): dimension of the graphical window
#Case with two coordinates
if(length(AP) == 2){
AP.list <- unique(spots.table$AP)
AP.list <- sort(AP.list[AP.list >= min(AP) & AP.list <= max(AP)], decreasing = TRUE)
for (i in AP.list)
plot.2d.atlas.expr(gene,spots.table,st.data,i,NULL,min.expr,cex,TRUE,text.plot,log.scale, normalized.scale)
}else{
if (!is.null(AP))
AP <- select.closest.AP(AP, spots.table)
#Case where plotting virtual slice (both dupplicates at a specific coordinates)
if (is.null(slice_index))
slice_index <- unique(spots.table[spots.table$AP == AP,'slice_index'])
old.slice.index <- slices.table[slice_index,'slice_old_id']
#Create a function to generate a continuous color palette
rbPal <- colorRampPalette(c('white','red'))
dataset <- get.expression.slice(gene, spots.table, st.data, min.expr, log.scale, normalized.scale, slice_index, rbPal)
AP <- dataset[[1]]$AP[1]
#Crashing case
if (length(old.slice.index) > 1){
for(j in 2:length(old.slice.index)){
if(dataset[[j]]$AP[1] != AP)
stop('Unxpected multiple AP coordinates')
}
}
m.val <- dataset[[(length(dataset) -1)]]
min.val <- dataset[[length(dataset)]]
#Selecting proper infos from the atlas
AP.infos <- atlas.spots$outlines$AP[atlas.spots$outlines$AP$AP == AP,]
outlines <- atlas.spots$outlines$outlines[[AP.infos$i]]
width <- AP.infos$xmax - AP.infos$xmin
height <- AP.infos$ymax - AP.infos$ymin
margin.x <- 0.4 * width
margin.y <- 0.05 * height
xMin <- AP.infos$xmin - margin.x
xMax <- AP.infos$xmax + margin.x
yMin <- AP.infos$ymin - margin.y
yMax <- AP.infos$ymax + margin.y
#Initializing the plot
if (new.window)
open.new.window(os.windows = os.windows, win.dim, X11.type = X11.type)
if (text.plot){
title.str <- paste("Bregma:", AP,"mm\nSlices: ",paste(slice_index, collapse = ' & '),'\nGene:',gene)
}else{
title.str <- ''
}
plot(c(xMin, xMax), c(yMin, yMax), ylim = c(yMax, yMin),
xlim = c(xMin, xMax), asp = 1, axes = F, xlab = "", ylab = "",
col = 0, main = title.str, font.main = 1)
if (text.plot){
mtext("Dorso-ventral (mm)", side = 2, line = 1.5)
mtext("Medio-lateral (mm)", side = 1, line = -1.5)
}
cnt.points <- 0
for (j in 1:length(outlines))
cnt.points <- cnt.points + length(outlines[[j]]$x)
df.points.plot <- matrix(ncol = 2, nrow = cnt.points+(length(outlines)-1))
curs <- 1
for (j in 1:length(outlines)){
o <- outlines[[j]]
df.points.plot[curs:(curs+length(o$x)-1),] <- cbind(as.numeric(o$x), as.numeric(o$y))
curs <- curs + length(o$x) + 1
#Colors for fiber and ventricle
if (as.character(o$col) == '#aaaaaa')
col = 'black'
else if (as.character(o$col) == '#cccccc')
col = 'darkgray'
else
col = NA
if(!is.na(col))
polygon(o$x,o$y,border = "black",col = col)
}
lines(df.points.plot)
#Plotting the spots
for(i in 1:length(old.slice.index)){
#Selecting the spots registered into the atlas
cur.soma <- atlas.spots$spots[rownames(dataset[[i]]),]
points(cur.soma$somaX, cur.soma$somaY, pch=21, bg = dataset[[i]]$col,cex = cex)
if (!is.null(spots.to.color))
points(cur.soma[spots.to.color,'somaX'], cur.soma[spots.to.color, 'somaY'], col = 'blue', pch=21, cex = cex, lwd = 2)
if (!is.null(spots.to.color.2))
points(cur.soma[spots.to.color.2,'somaX'], cur.soma[spots.to.color.2, 'somaY'], col = 'chartreuse3', pch=21, cex = cex, lwd = 2)
}
#Legend plotting
leg.ys <- (yMin+yMax)/2-0.25*height
leg.ye <- (yMin+yMax)/2+0.25*height
y.list <- seq(from = leg.ys, to = leg.ye, length.out = 101)
col.vect <- rev(rbPal(100))
x.list.mi <- numeric((length(y.list)-1))
x.list.ma <- numeric((length(y.list)-1))
y.list.mi <- numeric((length(y.list)-1))
y.list.ma <- numeric((length(y.list)-1))
for(i in 1:(length(y.list)-1)){
x.list.mi[[i]] <- xMax-0.25*width
y.list.mi[[i]] <- y.list[i]
x.list.ma[[i]] <- xMax-0.2*width
y.list.ma[[i]] <- y.list[i+1]
}
rect(x.list.mi,y.list.mi,x.list.ma,y.list.ma,col=col.vect, border = NA)
rect(xMax-0.25*width,leg.ys,xMax-0.2*width,leg.ye)
text(xMax-0.2*width,leg.ye, sprintf('%.0f',min.val), cex = cex,pos = 4)
text(xMax-0.2*width,leg.ys,sprintf('%.0f',m.val),cex = cex,pos = 4)
if(log.scale)
text(xMax-0.225*width,leg.ys,'Log scale',cex = cex,pos=3)
if(normalized.scale)
str.s <- 'Normalized scale'
else
str.s <- 'Number of reads'
text(xMax-0.15*width,(yMin + yMax)/2,str.s,cex = cex,pos=3,srt = -90)
if (!is.null(cluster.name.1))
text(xMin,yMin, cluster.name.1, pos=4, col = 'blue')
if (!is.null(cluster.name.2))
text(xMin,yMin + 0.025*height, cluster.name.2, pos=4, col = 'chartreuse3')
}
}
#Plotting expression levels on a slice
plot.2d.slice.expr <- function(gene,spots.table, st.data, slice_index, min.expr = 0, cex = 1, text.plot = TRUE, new.window = TRUE, log.scale = FALSE, normalized.scale = FALSE, win.dim = c(14,14), res.ds = 10, path.data.dir, os.windows = FALSE, X11.type = NULL){
#gene(str): gene to plot
#spots.table(df)
#st.data(df)
#slice_index(str): Section to plot
#min.expr(float): Expression threhsold to plot spot. If between 0 and 1, defines a percentile. Otherwise, defines a number of reads (not log normalized)
#cex(float): dimension of text
#new.window(bool): force new window
#text.plot(bool): show the legend and title
#log.scale(bool): use a logarithmic scale instead
#normalized.scale(bool): all expression are a value between 0 and 100 (maximum expression across the whole brain)
#win.dim(numeric vector): dimension of the graphical window
#res.ds(integer): downsample factor for plotting the HE picture faster
#Create a function to generate a continuous color palette
rbPal <- colorRampPalette(c('white','red'))
dataset <- get.expression.slice(gene, spots.table, st.data, min.expr, log.scale, normalized.scale, slice_index, rbPal)
spots.cur <- dataset[[1]]
spots.cur$radius <- spots.table[rownames(spots.cur),'radius']
m.val <- dataset[[2]]
min.val <- dataset[[3]]
old.slice.index <- slices.table[slice_index,'slice_old_id']
rotation <- slices.table[slice_index,'rotation']
HE.path <- sprintf('%s/pictures/HE_%s.jpg',path.data.dir,old.slice.index)
HE.original <- readJPEG(HE.path)
HE.resolution.x <- dim(HE.original)[2]
HE.resolution.y <- dim(HE.original)[1]
HE.downscale <- HE.original[seq(1,HE.resolution.y,by=res.ds),seq(1,HE.resolution.x,by=res.ds),]
if(new.window)
open.new.window(os.windows = os.windows, win.dim, X11.type = X11.type)
if(rotation == 90){
xlim <- c(0,HE.resolution.x)
ylim <- c(2*HE.resolution.y,HE.resolution.y)
angle <- -rotation
spots.cur$y <- spots.cur$y + HE.resolution.y
}
if(rotation == -90){
xlim <- c(-HE.resolution.x,0)
ylim <- c(HE.resolution.y,0)
angle <- -rotation
spots.cur$x <- spots.cur$x - HE.resolution.y
spots.cur$y <- spots.cur$y + (abs(HE.resolution.x - HE.resolution.y))
}
if(rotation == 0){
xlim <- c(0,HE.resolution.x)
ylim <- c(HE.resolution.y,0)
angle <- 0
}
if (text.plot)
title.str <- paste("Bregma:", spots.cur$AP[1],"mm\nSlices: ",paste(slice_index, collapse = ' & '),'\nGene:',gene)
else
title.str <- ''
width <- abs(xlim[1] - xlim[2])
plot(c(xlim[1],xlim[2] + 0.2*width), ylim, xlim = c(xlim[1],xlim[2] + 0.4*width), ylim = ylim,
asp = 1, axes = F, xlab = "", ylab = "",
col = 0, main = title.str, font.main = 1)
rasterImage(as.raster(HE.downscale), xleft = 1, xright = HE.resolution.x, ybottom = HE.resolution.y, ytop = 1, angle=angle, asp = 1)
symbols(spots.cur$x,spots.cur$y,circles = spots.cur$radius, bg = spots.cur$col, inches=FALSE, asp = 1,fg='black',add = TRUE, lwd = 0.2, asp = 1)
leg.ys <- mean(ylim)-(0.2*abs(ylim[1]-ylim[2]))
leg.ye <- mean(ylim)+(0.2*abs(ylim[1]-ylim[2]))
y.list <- seq(from = leg.ys, to = leg.ye, length.out = 101)
col.vect <- rev(rbPal(100))
xMax <- max(xlim)
x.list.mi <- numeric((length(y.list)-1))
x.list.ma <- numeric((length(y.list)-1))
y.list.mi <- numeric((length(y.list)-1))
y.list.ma <- numeric((length(y.list)-1))
for(i in 1:(length(y.list)-1)){
x.list.mi[[i]] <- xMax+0.1*width
y.list.mi[[i]] <- y.list[i]
x.list.ma[[i]] <- xMax+0.15*width
y.list.ma[[i]] <- y.list[i+1]
}
rect(x.list.mi,y.list.mi,x.list.ma,y.list.ma,col=col.vect, border = NA)
rect(xMax+0.1*width,leg.ys,xMax+0.15*width,leg.ye)
text(xMax+0.15*width,leg.ye,sprintf('%.0f',min.val),cex = cex,pos = 4)
text(xMax+0.15*width,leg.ys,sprintf('%.0f',m.val),cex = cex,pos = 4)
if(log.scale)
text(xMax+0.125*width,leg.ys,'Log scale',cex = cex,pos=3)
if(normalized.scale)
str.s <- 'Normalized scale'
else
str.s <- 'Number of reads'
text(xMax+0.17*width,mean(ylim),str.s,cex = cex,pos=3,srt = -90)
}
open.3d.plot <- function(p){
path <- paste(tempfile(),'html',sep='.')
p <- htmlwidgets::saveWidget(p, path, selfcontained = FALSE)
browseURL(path)
}
## ----------- Window manager -----------
#Open a new window depending on the OS
open.new.window <- function(os.windows, win.dim, X11.type){
if (os.windows){
win.graph(width = win.dim[1], height = win.dim[2])
}else{
if (is.null(X11.type))
quartz(width = win.dim[1], height = win.dim[2])
else
X11(width = win.dim[1], height = win.dim[2], type=X11.type)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.