R/plotTrackOverlay_Dcoef.R

In sojourner: Statistical analysis of single molecule trajectories

#### 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(folder=folder, input=3)
##' 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 = TRUE)
    }
    
    ## 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(seq_along(trackll))) {
        
        x = (2:5) * t.interval
        dimension = 2
        trackl = list()
        
        for (j in c(seq_along(trackll[[i]]))) {
            track = trackll[[i]][[j]]
            msd.n = msd.track(track, dt = dt, at.dt = FALSE)
            
            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 = vapply(trackl, function(x) {
            x[[2]]
        },double(1))
        Dcoef[!is.finite(Dcoef)] <- NA
        trackl = trackl[!is.na(Dcoef)]
        Dcoef = vapply(trackl, function(x) {
            x[[2]]
        }, double(1))
        trackl = trackl[which(Dcoef >= (Dcoef.range[1]))]
        Dcoef = vapply(trackl, function(x) {
            x[[2]]
        }, double(1))
        trackl = trackl[which(Dcoef <= (Dcoef.range[2]))]
        Rsquare = vapply(trackl, function(x) {
            x[[3]]
        },double(1))
        trackl = trackl[which(vapply(Rsquare, function(x) {
            x >= rsquare
        }, logical(1)))]

        Dcoef = vapply(trackl, function(x) {
            x[[2]]
        },double(1))

        ## 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(seq_along(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 seq_len(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 = 0.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()
    
}