Merge pull request #1 from majpark21/master

Power Spectral Density

global.R

@@ -6,6 +6,7 @@ source('modules/dispStats.R')
 source('modules/dispTrackStats.R')
 source('modules/trajPlot.R')
 source('modules/trajRibbonPlot.R')
+source('modules/trajPsdPlot.R')
 source('modules/boxPlot.R')
 source('modules/tabAUC.R')
 source('modules/clDistPlot.R')

modules/auxfunc.R

@@ -53,6 +53,7 @@ FCSVOUTLIERS = 'outliers.csv'
 FCSVTCCLEAN  = 'tCoursesSelected_clean.csv'
 FPDFTCMEAN   = "tCoursesMeans.pdf"
 FPDFTCSINGLE = "tCourses.pdf"
+FPDFTCPSD    = 'tCoursesPsd.pdf'
 FPDFBOXAUC   = 'boxplotAUC.pdf'
 FPDFBOXTP    = 'boxplotTP.pdf'
 FPDFSCATTER  = 'scatter.pdf'
@@ -206,6 +207,57 @@ LOCcalcTrajCI = function(in.dt, in.col.meas, in.col.by = NULL, in.type = c('norm
     return(loc.dt)
 }
 
+
+#' Calculate the power spectrum density for time-series
+#'
+#' @param in.dt Data table in long format
+#' @param in.col.meas Name of the column with the measurement
+#' @param in.col.id Name of the column with the unique series identifier
+#' @param in.col.by Column names for grouping (default NULL - no grouping). PSD of individual trajectories will be averaged within a group.
+#' @param in.method Name of the method for PSD estimation, must be one of c("pgram", "ar"). Default to "pgram*.
+#' @param in.return.period Wheter to return densities though periods (1/frequencies) instead of frequencies.
+#' @param ... Other paramters to pass to stats::spectrum()
+#'
+#' @return Datatable with columns: (frequency or period), spec (the density) and grouping column
+#' @export
+#' @import data.table
+#'
+#' @examples
+LOCcalcPSD <- function(in.dt,
+                    in.col.meas,
+                    in.col.id,
+                    in.col.by,
+                    in.method = "pgram",
+                    in.return.period = TRUE,
+                    ...){
+  require(data.table)
+  # Method "ar" returns $spec as matrix whereas "pgram" returns a vector, custom function to homogenze output format
+  mySpectrum <- function(x, ...){
+    args_spec <- list(x=x, plot=FALSE)
+    inargs <- list(...)
+    args_spec[names(inargs)] <- inargs
+    out <- do.call(spectrum, args_spec)
+    out$spec <- as.vector(out$spec)
+    return(out)
+  }
+  if(!in.method %in% c("pgram", "ar")){
+    stop('Method should be one of: c("pgram", "ar"')
+  }
+  dt_spec <-  in.dt[, (mySpectrum(get(in.col.meas), plot = FALSE, method = in.method)[c("freq", "spec")]), by = in.col.id]
+  dt_group <- in.dt[, .SD[1, get(in.col.by)], by = in.col.id]
+  setnames(dt_group, "V1", in.col.by)
+  dt_spec <- merge(dt_spec, dt_group, by = in.col.id)
+  dt_agg <- dt_spec[, .(spec = mean(spec)), by = c(in.col.by, "freq")]
+  if(in.return.period){
+    dt_agg[, period := 1/freq]
+    dt_agg[, frequency := NULL]
+  } else {
+    setnames(dt_agg, "freq", "frequency")
+  }
+  return(dt_agg)
+}
+
+
 #' Generate synthetic CellProfiler output with single cell time series
 #'
 #'
@@ -644,7 +696,34 @@ LOCplotTrajRibbon = function(dt.arg, # input data table
   return(p.tmp)
 }
 
-
+# Plot average power spectrum density per facet
+LOCplotPSD <- function(dt.arg, # input data table
+                    x.arg, # string with column name for x-axis
+                    y.arg, # string with column name for y-axis
+                    group.arg=NULL, # string with column name for grouping time series (here, it's a column corresponding to grouping by condition)
+                    xlab.arg = x.arg,
+                    ylab.arg = y.arg,
+                    col.arg = NULL){
+  require(ggplot2)
+  if(length(setdiff(c(x.arg, y.arg, group.arg), colnames(dt.arg))) > 0){
+    stop(paste("Missing columns in dt.arg: ", setdiff(c(x.arg, y.arg, group.arg), colnames(dt.arg))))
+  }
+  p.tmp <- ggplot(dt.arg, aes_string(x=x.arg, y=y.arg)) +
+    geom_line() +
+    geom_rug(sides="b", alpha = 1, color = "lightblue") +
+    facet_wrap(group.arg) +
+    labs(x = xlab.arg, y = ylab.arg)
+  
+  if (is.null(col.arg)) {
+    p.tmp = p.tmp +
+      scale_color_discrete(name = '')
+  } else {
+    p.tmp = p.tmp +
+      scale_colour_manual(values = col.arg, name = '')
+  }
+  
+  return(p.tmp)
+}
 
 # Plots a scatter plot with marginal histograms
 # Points are connected by a line (grouping by cellID)

modules/tabClHier.R

@@ -172,6 +172,9 @@ clustHierUI <- function(id, label = "Hierarchical CLustering") {
       tabPanel('Time-courses',
                modTrajPlotUI(ns('modPlotHierTraj'))),
       
+      tabPanel('PSD',
+               modPSDPlotUI(ns('modPlotHierPsd'))),
+      
       tabPanel('Cluster dist.',
                modClDistPlotUI(ns('hierClDistPlot'), 'xxx'))
       
@@ -528,13 +531,22 @@ clustHier <- function(input, output, session, in.data4clust, in.data4trajPlot, i
            '.pdf')
   })
   
+  createFnamePsdPlot = reactive({
+    
+    paste0('clust_hierch_tCoursesPsd_',
+           s.cl.diss[as.numeric(input$selectPlotHierDiss)],
+           '_',
+           s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], 
+           '.pdf')
+  })
+  
   createFnameDistPlot = reactive({
     
     paste0('clust_hierch_clDist_',
            s.cl.diss[as.numeric(input$selectPlotHierDiss)],
            '_',
-           s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.pdf')  })
-  
+           s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.pdf')  
+  })
   
   
   #  Hierarchical - Heat Map - download pdf
@@ -556,6 +568,13 @@ clustHier <- function(input, output, session, in.data4clust, in.data4trajPlot, i
              in.facet.color = getClColHier,
              in.fname = createFnameRibbonPlot)
   
+  # plot cluster PSD
+  callModule(modPSDPlot, 'modPlotHierPsd',
+             in.data = data4trajPlotCl,
+             in.facet = 'cl',
+             in.facet.color = getClColHier,
+             in.fname = createFnamePsdPlot)
+
   # plot distribution barplot
   callModule(modClDistPlot, 'hierClDistPlot', 
              in.data = data4clDistPlot,

modules/trajPsdPlot.R

+require(DT)
+
+ # UI ----
+modPSDPlotUI =  function(id, label = "Plot PSD of average trajectory.") {
+  ns <- NS(id)
+  
+  tagList(
+    fluidRow(
+      column(
+        3,
+        checkboxInput(ns('chBplotTrajInt'), 'Interactive Plot'),
+        actionButton(ns('butPlotTraj'), 'Plot!')
+      ),
+      column(
+        3,
+        selectInput(ns('inPSDxchoice'), 'Xaxis:', list('Period'= TRUE, 'Frequency'= FALSE)),
+        radioButtons(ns('rBPSDmethod'), 'Method for PSD estimation:', list('Smoothed Fourier' = 'pgram', 'AR Fit' = 'ar'))
+      ),
+      column(
+        3,
+        selectInput(ns('inPSDlogtype'), 'Log function:', list('log2'= 'log2', 'log10'= 'log10', 'ln'= 'log')),
+        checkboxGroupInput(ns('chBGPSDlogaxis'), 'Log the axis:', list('x' = 'x', 'y' = 'y'), inline = TRUE)
+      ),
+      column(
+        3,
+        numericInput(
+          ns('inPlotTrajWidth'),
+          'Width [%]:',
+          value = 100,
+          min = 10,
+          max = 100,
+          width = '100px',
+          step = 10
+        ),
+        numericInput(
+          ns('inPlotTrajHeight'),
+          'Height [px]:',
+          value = 800,
+          min = 100,
+          width = '100px',
+          step = 50
+        )
+      )
+    ),
+    uiOutput(ns('uiPlotTraj')),
+    br(),
+    modTrackStatsUI(ns('dispTrackStats')),
+    
+    downPlotUI(ns('downPlotTraj'), "Download PDF")
+  )
+}
+
+# Server ----
+
+modPSDPlot = function(input, output, session, 
+                             in.data,
+                             in.fname,
+                             in.facet = 'group', 
+                             in.facet.color = NULL
+                             ) {
+  
+  ns <- session$ns
+  
+  output$uiPlotTraj = renderUI({
+    if (input$chBplotTrajInt)
+      plotlyOutput(
+        ns("outPlotTrajInt"),
+        width = paste0(input$inPlotTrajWidth, '%'),
+        height = paste0(input$inPlotTrajHeight, 'px')
+      ) else
+        plotOutput(
+          ns("outPlotTraj"),
+          width = paste0(input$inPlotTrajWidth, '%'),
+          height = paste0(input$inPlotTrajHeight, 'px')
+        )
+  })
+  
+  
+  callModule(modTrackStats, 'dispTrackStats',
+             in.data = in.data)
+  
+  
+  output$outPlotTraj <- renderPlot({
+    
+    loc.p = plotTraj()
+    if(is.null(loc.p))
+      return(NULL)
+    
+    return(loc.p)
+  })
+  
+  
+  output$outPlotTrajInt <- renderPlotly({
+    # This is required to avoid
+    # "Warning: Error in <Anonymous>: cannot open file 'Rplots.pdf'"
+    # When running on a server. Based on:
+    # https://github.com/ropensci/plotly/issues/494
+    if (names(dev.cur()) != "null device")
+      dev.off()
+    pdf(NULL)
+    
+    loc.p = plotTraj()
+    if(is.null(loc.p))
+      return(NULL)
+    
+    return(plotly_build(loc.p))
+  })
+  
+  
+  
+  # PSD plot - download pdf
+  callModule(downPlot, "downPlotTraj", 
+             in.fname = in.fname,
+             plotTraj, TRUE)
+  
+  plotTraj <- function() {
+    cat(file = stderr(), 'plotPSD: in\n')
+    locBut = input$butPlotTraj
+    
+    if (locBut == 0) {
+      cat(file = stderr(), 'plotPSD: Go button not pressed\n')
+      
+      return(NULL)
+    }
+    
+    # check if main data exists
+    loc.dt = isolate(in.data())
+    
+    cat("plotPSD: on to plot\n\n")
+    if (is.null(loc.dt)) {
+      cat(file = stderr(), 'plotPSD: dt is NULL\n')
+      return(NULL)
+    }
+    
+    cat(file = stderr(), 'plotPSD: dt not NULL\n')
+    
+    
+    
+    # Future: change such that a column with colouring status is chosen by the user
+    # colour trajectories, if dataset contains mid.in column
+    # with filtering status of trajectory
+    if (sum(names(loc.dt) %in% 'mid.in') > 0)
+      loc.line.col.arg = 'mid.in'
+    else
+      loc.line.col.arg = NULL
+    
+    # select every other point for plotting (fixed for PSD because lead to false interpretation of PSD)
+    loc.dt = loc.dt[, .SD[seq(1, .N, 1)], by = id]
+
+    # check if columns with XY positions are present
+    if (sum(names(loc.dt) %like% 'pos') == 2)
+      locPos = TRUE
+    else
+      locPos = FALSE
+    
+    # check if column with ObjectNumber is present
+    if (sum(names(loc.dt) %like% 'obj.num') == 1)
+      locObjNum = TRUE
+    else
+      locObjNum = FALSE
+    
+    
+    # If in.facet.color present,
+    # make sure to include the same number of colours in the palette,
+    # as the number of groups in dt.
+    # in.facet.color is typically used when plotting time series within clusters.
+    # Then, the number of colours in the palette has to be equal to the number of clusters (facetted according to in.facet variable).
+    # This might differ if the user selects manually clusters to display.
+    if (is.null(in.facet.color)) 
+      loc.facet.col = NULL 
+    else {
+      # get group numbers in dt; 
+      # loc.dt[, c(in.facet), with = FALSE] returns a data table with a single column
+      # [[1]] at the end extracts the first column and returns as a vector
+      loc.groups = unique(loc.dt[, c(in.facet), with = FALSE][[1]])
+      
+      # get colour palette
+      # the length is equal to the number of groups in the original dt.
+      # When plotting time series within clusters, the length equals the number of clusters.
+      loc.facet.col = in.facet.color()$cl.col
+      loc.facet.col = loc.facet.col[loc.groups]
+    }
+    
+
+    loc.dt.aggr <- LOCcalcPSD(in.dt = loc.dt,
+                              in.col.meas = 'y',
+                              in.col.id = 'id',
+                              in.col.by = in.facet,
+                              in.method = input$rBPSDmethod,
+                              in.return.period = input$inPSDxchoice
+                              )
+    loc.dt.aggr[, (in.facet) := as.factor(get(in.facet))]
+    
+    x_arg <- ifelse('period' %in% colnames(loc.dt.aggr), 'period', 'frequency')
+    x_arg_str <- paste0(toupper(substr(x_arg, 1, 1)), tolower(substring(x_arg, 2)))  # capitalized
+    p.out <- LOCplotPSD(dt.arg = loc.dt.aggr,
+                        x.arg = x_arg,
+                        y.arg = 'spec',
+                        group.arg = in.facet,
+                        col.arg = loc.facet.col,
+                        xlab.arg = x_arg_str,
+                        ylab.arg = '') +
+      LOCggplotTheme(in.font.base = PLOTFONTBASE, 
+                     in.font.axis.text = PLOTFONTAXISTEXT, 
+                     in.font.axis.title = PLOTFONTAXISTITLE, 
+                     in.font.strip = PLOTFONTFACETSTRIP, 
+                     in.font.legend = PLOTFONTLEGEND)
+    if("x" %in% input$chBGPSDlogaxis){
+      p.out <- p.out + scale_x_continuous(trans = input$inPSDlogtype)
+    }
+    if("y" %in% input$chBGPSDlogaxis){
+      p.out <- p.out + scale_y_continuous(trans = input$inPSDlogtype)
+    }
+    return(p.out)
+  }
+}
\ No newline at end of file

server.R

@@ -823,6 +823,11 @@ shinyServer(function(input, output, session) {
              in.data.stim = data4stimPlot,
              in.fname = function() {return(FPDFTCSINGLE)})
   
+  # Trajectory plotting - PSD
+  callModule(modPSDPlot, 'modPSDPlot',
+             in.data = data4trajPlotNoOut,
+             in.fname = function() {return(FPDFTCPSD)})
+  
   
   # Tabs ----
   ###### AUC calculation and plotting

ui.R

@@ -118,6 +118,12 @@ shinyUI(fluidPage(
               uiOutput('varSelHighlight'),
               br(),
               modTrajPlotUI('modTrajPlot')
+            ),
+            
+            tabPanel(
+              "Power Spectral Density",
+              br(),
+              modPSDPlotUI('modPSDPlot')
             )
           )
         ),