R/plotFunctions.R

In cantin-ortiz/stExpressionViewer: Plotting gene expression from the ST Data.

## ----------- 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)
  }
}