sojourner: Single Molecule Tracking

#### plotTrackOverlay.Dcoef.R
#### Wu Lab, Johns Hopkins University
#### Author: Xiaona Tang
#### Date: Dec 20, 2017

## plotTrackOverlay.Dcoef-methods
##
###############################################################################
##' @name plotTrackOverlay.Dcoef
##' @aliases plotTrackOverlay.Dcoef
##' @title plotTrackOverlay.Dcoef
##' @rdname plotTrackOverlay.Dcoef-methods
##' @docType methods
##' @description Plot track/trajectory overlays from a list of files in a folder with color coded 
##'              by the Diffusion Coefficient (Dcoef) of each track/trajectory. 
##'              
##' @usage
##' plotTrackOverlay.Dcoef(trackll=trackll,scale=128,dt=6,filter=c(min=7,max=Inf),
##' resolution=0.107,rsquare=0.8,t.interval=0.01,Dcoef.range=c(-6,2),
##' color=c("blue", "white", "red"),folder=NULL,file.No=0,line.width=0.1)
##'
##' @param trackll trajectory list generated by createTrackll() and processing. if NULL, user will be prompted to enter the trackll name.
##' @param scale The pixel scale of image data.
##' @param dt Time inbervals.
##' @param filter Filter the tracks by step/frame number (length). Only tracks pass through filter will be selected.
##' @param resolution ratio of pixel to uM.
##' @param rsquare R square filter on Dcoef results. Default to be 0.8. Set value to 0 if rsquare filter is not desired.
##' @param t.interval Time interval for image aquisition. Default 0.01 sec (10ms).
##' @param Dcoef.range Select tracks whose Dcoef is within the range. Other tracks will not be plotted.
##' @param color colors to interpolate to build the color gradient/ramp.
##' @param folder Path to provide nuclear glow backround. Images should end with "_Nuclei.tif". Default is NULL and not adding background.
##' @param file.No Select file(s) in the folder to plot. Default 0 for plotting all files in the folder.
##' @param line.width Line width for plotted tracks, default is 0.1.
##' @return
##' \itemize{
##' \item{PDF:} One PDF file with one plot on each page.
##' }
##' @details Plot track/trajectory overlays from a list of files in a folder with color coded 
##'          by the Diffusion Coefficient (Dcoef) of each track/trajectory. The range of Dcoef can be defined. Nuclei backgrounds can be added to the plots.
##'          Tracks from each file will be plot in one plot, i.e., one plot per file. Colors are flexible.
##'              
##'          Upon running of the function, users will be prompted to input the name of the track list (trackll).
##'          Input un-merged trackll and the plotting will start. 
##'          The Dcoef is calculated by "static" method, which stabilizes the number of time lags used for fitting using time lag 2~ 5 despite the total time lags measured.
##'
##'
##' @examples
##'
##' # Generate trackll, and process, 
##' # e.g. mask region of interest, tracks from multiple files should not be merged.
##' folder=system.file("extdata","HSF",package="sojourner")
##' trackll=createTrackll(interact=FALSE,folder,input=2, cores = 2)
##' trackll=maskTracks(folder,trackll)
##' 
##' # Plot track overlays,
##' plotTrackOverlay.Dcoef(trackll=trackll,scale=128,dt=6,filter=c(min=7,max=Inf),
##' resolution=0.107,rsquare=0.8,t.interval=0.01,Dcoef.range=c(-6,2),
##' color=c("blue", "white", "red"),folder=NULL,file.No=0,line.width=0.1)
##' plotTrackOverlay.Dcoef(trackll=trackll,color=c("red", "yellow"),folder=folder,file.No=0)
##' plotTrackOverlay.Dcoef(trackll=trackll,color=c("red", "yellow"),folder=folder,file.No=c(1,2))

##' @export plotTrackOverlay.Dcoef

#####################################################################################
#####################################################################################


## Function for plotting trajectory overlay with color coded by Dcoef
.plotTrackOverlay.Dcoef<-function(trackll=trackll,scale=128,dt=6,filter=c(min=7,max=Inf),resolution=0.107,rsquare=0.8,
                          t.interval=0.01,Dcoef.range=c(-6,2),color=c("blue", "white", "red"),
                          folder=NULL,file.No=0,line.width=0.1){

  ## Import trackll (un-merged) information
  if (length(file.No)>1&min(file.No) == 0){
    stop("Wrong file.No input. Ceased.",call. = FALSE, domain = NULL)
  }
  
  if(is.null(trackll)){
    trackll.label<-readline(cat("Enter the un-merged trackll you want to plot:    "))
    trackll<-get(paste(trackll.label))
  }
  
  if (length(file.No)>=1&file.No[1]>0){
    trackll<-trackll[file.No]
  }
  
  trackll<-filterTrack(trackll,filter=c(min=dt+1,max=Inf))
  trackll<-filterTrack(trackll,filter=filter)
  
  ## Set plot area as black background and white frontground.
  oldpar <- par
  par(mar=c(3, 4, 4, 4),xpd=FALSE)
  
  par(mfrow=c(1,1),bg="black",fg="white")
  
  cat("Plotting...A PDF file will be output in the working directory.\n")
  
  ## Read in nuclei background data if the image folder is provided.
  if(!is.null(folder)){
    nuclei.lst=list.files(path=folder,pattern="_Nuclei.tif",full.names=T)
  }
  
  ## Generate color gradient/ramp for trajectory plotting.
  #cl=grDevices::heat.colors(ceiling(max(Dcoef)-min(Dcoef))*10)
  #cl=rev(cl)
  #cl=colorspace::diverge_hsv(max(dense))
  #cl=colorspace::diverge_hcl(max(dense))
  cl <- colorRampPalette(color)(n = 101)
  
  #Get trackl information in the un-merged trackll and filter by Dcoef range and rsquare.
  for (i in c(1:length(trackll))){

    x=(2:5)*t.interval
    dimension=2
    trackl=list()
    
    for (j in c(1:length(trackll[[i]]))){
      track=trackll[[i]][[j]]
      msd.n=msd.track(track, dt=dt, at.dt=F)
      
      fit=lm(msd.n[2:5]~x)
      MSDslope=coefficients(fit)[2]/(2*dimension)
      Log.D.coef=log(MSDslope)
      MSDcorr=summary(fit)$r.squared
      #if((!is.na(Log.D.coef))&(MSDcorr>=rsquare))
      trackl[[j]]=list(track,Log.D.coef,MSDcorr)
    }
    #trackl=trackl[-which(sapply(trackl, is.null))]
    
    Dcoef=sapply(trackl,function(x){x[[2]]})
    Dcoef[!is.finite(Dcoef)] <- NA
    trackl=trackl[!is.na(Dcoef)]
    Dcoef=sapply(trackl,function(x){x[[2]]})
    trackl=trackl[which(Dcoef>=(Dcoef.range[1]))]
    Dcoef=sapply(trackl,function(x){x[[2]]})
    trackl=trackl[which(Dcoef<=(Dcoef.range[2]))]
    Rsquare=sapply(trackl,function(x){x[[3]]})
    trackl=trackl[which(sapply(Rsquare, function(x){x>=rsquare}))]
    
    Dcoef=sapply(trackl,function(x){x[[2]]})
    
    ## Plot trackoverlay for each file (trackl) in the un-merged trackll.
    plot.new()
    ## If image folder is provided, add nuclei background to the plot.
    if(!is.null(folder)){
      #library(EBImage)
      img=EBImage::readImage(nuclei.lst[[i]])
      d=img@.Data
      raster.img <- as.raster(t(d[,,1]))
      plot.window(xlim=c(0,scale),ylim=c(0,scale),xaxs = "i", yaxs = "i")
      plot(raster.img,add = TRUE)
    }
    mtext(gsub(".mat","",names(trackll[i])),side=3,line=0.5,cex=2,col.main="white") 
    
    ## Rescale the plotting area with resolution, changing the units from pixel to um.
    plot.window(xlim=c(0,scale*resolution),ylim=c(0,scale*resolution),xaxs = "i", yaxs = "i")
    axis(1,cex.axis=1,col.axis="white")
    axis(2,cex.axis=1,col.axis="white")
    mtext(expression(paste("X (",mu,"m)")),side=1,line=2,cex.lab=1,col="white")
    mtext(expression(paste("Y (",mu,"m)")),side=2,line=2,cex.lab=1,col="white")
    box()
    for(k in c(1:length(trackl))){
      lines(trackl[[k]][[1]]$x*resolution,(128-trackl[[k]][[1]]$y)*resolution,
            col=cl[(Dcoef[k]-Dcoef.range[1])/(Dcoef.range[2]-Dcoef.range[1])*100+1],lwd=line.width)
      
    }
    
    ## Add color gradient legend to the right edge of each plot.
    .legend.col(col = cl, lev = Dcoef.range)
    
    ## Add parameters and track number (n) as text legend to the topright corner of each plot.
    legend=c(paste("dt = ",dt),as.expression(bquote('Rsquare'>=.(rsquare))), 
             paste("n = ",length(Dcoef)),paste('Dcoef range [',Dcoef.range[1],",",Dcoef.range[2],"]"))
    
    temp <- legend("topright", legend = c(" ", " "," "," "),
                   text.width = strwidth("Rsquare>=0.8"),
                   xjust = 1, yjust = 1,bty="n",cex=2)
    text(temp$rect$left + temp$rect$w, temp$text$y,
         legend,col="white",pos = 2,cex=2)
    
  }
  
  ## Reset plotting area parameters.
  par(oldpar)
  par(mfrow=c(1,1),bg="white",fg="black")
}


## Function for adding color gradient legend.
## Refernce: by Aur?lien. 
## From: https://aurelienmadouasse.wordpress.com/2012/01/13/legend-for-a-continuous-color-scale-in-r/
.legend.col <- function(col, lev){
  
  opar <- par
  
  n <- length(col)
  
  bx <- par("usr")
  
  box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000,
              bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
  box.cy <- c(bx[3], bx[3])
  box.sy <- (bx[4] - bx[3]) / n
  
  xx <- rep(box.cx, each = 2)
  
  par(xpd = TRUE)
  for(i in 1:n){
    
    yy <- c(box.cy[1] + (box.sy * (i - 1)),
            box.cy[1] + (box.sy * (i)),
            box.cy[1] + (box.sy * (i)),
            box.cy[1] + (box.sy * (i - 1)))
    polygon(xx, yy, col = col[i], border = col[i])
    
  }
  par(new = TRUE)
  plot(0, 0, type = "n",
       ylim = c(min(lev), max(lev)),
       yaxt = "n", ylab = "",
       xaxt = "n", xlab = "",
       frame.plot = FALSE)
  axis(side = 4, las = 2, tick = FALSE, line = .25,col.axis="white")
  par (opar)
}



## Function for outputting the plots into one multipage PDF file in the working directory.
plotTrackOverlay.Dcoef<-function(trackll=trackll,scale=128,dt=6,filter=c(min=7,max=Inf),resolution=0.107,rsquare=0.8,
                         t.interval=0.01,Dcoef.range=c(-6,2),color=c("blue", "white", "red"),
                         folder=NULL,file.No=0,line.width=0.1){
  ## Output the plots into one PDF file in the working directory.
  pdf(paste("plotTrackOverlay.Dcoef.Heatmap--",format(Sys.time(),"%Y%m%d.%H%M%S"),".pdf",sep=""),width=11.7,height=11.7)

  .plotTrackOverlay.Dcoef(trackll=trackll,scale=scale,dt=dt,filter=filter,resolution=resolution,rsquare=rsquare,
                  t.interval=t.interval,Dcoef.range=Dcoef.range,color=color, 
                  folder=folder,file.No=file.No,line.width = line.width)
  
  dev.off()
  
}