Moved several functions to tca package

README.md

@@ -23,6 +23,20 @@ Following packages need to be installed in order to run the app locally:
 
 Install packages using `install.packages('name_of_the_package_from_the_list_above')` command in RStudio command line.
 
+**Additionally**, a time series analysis package need to be installed from [GitHub](https://github.com/dmattek/tca-package):
+
+```
+install_github("dmattek/tca-package")
+```
+
+The `install_github` function is available in *devtools* package:
+
+```
+install.packages("devtools")
+library(devtools)
+```
+
+
 #### Input file
 The app recognizes CSV (comma-separated values) files: data columns separated by a comma, floating point numbers using a dot (full-stop).
 

modules/auxfunc.R

@@ -95,7 +95,7 @@ help.text = c(
   Say, the main data file contains columns Metadata_Site and TrackLabel. 
   These two columns should be then selected in UI to form a unique cell ID, e.g. 001_0001 where former part corresponds to Metadata_Site and the latter to TrackLabel.',
   'Plotting and data processing requires a unique cell ID across entire dataset. A typical dataset from CellProfiler assigns unique cell ID (TrackLabel) within each field of view (Metadata_Site).
-                   Therefore, a unique ID is created by concatenating these two columns. If the dataset already contains a unique ID, check this box and select a single column only.'
+                   Therefore, a unique ID is created by concatenating these two columns. If the dataset already contains a unique ID, UNcheck this box and select a single column only.'
 )
 
 
@@ -317,41 +317,6 @@ myGgplotTraj = function(dt.arg, # data table
 }
 
 
-userDataGen <- function() {  
-  cat(file=stderr(), 'userDataGen: in\n')
-  
-  locNtp = 60
-  locNtracks = 10
-  locNsites = 6
-  locNwells = 1
-  
-  x.rand.1 = c(rnorm(locNtp * locNtracks * locNsites * 1/3, 0.5, 0.1), rnorm(locNtp * locNtracks * locNsites * 1/3,   1, 0.2), rnorm(locNtp * locNtracks * locNsites * 1/3,  2, 0.5))
-  x.rand.2 = c(rnorm(locNtp * locNtracks * locNsites * 1/3, 0.25, 0.1), rnorm(locNtp * locNtracks * locNsites * 1/3, 0.5, 0.2),  rnorm(locNtp * locNtracks * locNsites * 1/3, 1, 0.2))
-
-  # add NA's for testing
-  x.rand.1[c(10,20,30)] = NA
-  
-  #  x.rand.3 = rep(rnorm(locNtracks, 2, 0.5), 1, each = locNtp)
-#  x.rand.4 = rep(rnorm(locNtracks, 1, 0.1), 1, each = locNtp)
-  
-#  x.arg = rep(seq(0, locNtp-1) / locNtp * 4 * pi, locNtracks * locNsites)
-  x.arg = rep(seq(1, locNtp), locNtracks * locNsites)
-  
-  dt.nuc = data.table(Metadata_Site = rep(1:locNsites, each = locNtp * locNtracks),
-                      Metadata_Well = rep(1:locNwells, each = locNtp * locNsites * locNtracks / locNwells),
-                      Metadata_RealTime = x.arg,
-                      objCyto_Intensity_MeanIntensity_imErkCor = x.rand.1,
-                      objNuc_Intensity_MeanIntensity_imErkCor  = x.rand.2,
-                      objNuc_Location_X = runif(locNtp * locNtracks * locNsites, min = 0, max = 1),
-                      objNuc_Location_Y = runif(locNtp * locNtracks * locNsites, min = 0, max = 1),
-#                      objCyto_Intensity_MeanIntensity_imErkCor = x.rand.3 + ifelse(x.arg < 4, 0, 1) / x.rand.3,
-#                      objNuc_Intensity_MeanIntensity_imErkCor  = c(rnorm(locNtp * locNtracks * locNsites * 0.5, .25, 0.1), rnorm(locNtp * locNtracks * locNsites * 0.5, .5, 0.2)),
-                      TrackLabel = rep(1:(locNtracks*locNsites), each = locNtp))
-  
-  return(dt.nuc)
-}
-
-
 # Fast DTW computation
 fastDTW <-function (x)
 {
@@ -359,65 +324,6 @@ fastDTW <-function (x)
 }
 
 
-# Returns original dt with an additional column with normalized quantity.
-# The column to be normalised is given by 'in.meas.col'.
-# The name of additional column is the same as in.meas.col but with ".norm" suffix added.
-# Normalisation is based on part of the trajectory;
-# this is defined by in.rt.min and max, and the column with time in.rt.col.
-# Additional parameters:
-# in.by.cols - character vector with 'by' columns to calculate normalisation per group
-#              if NULL, no grouping is done
-# in.robust - whether robust measures should be used (median instead of mean, mad instead of sd)
-# in.type - type of normalization: z.score or mean (fi.e. old change w.r.t. mean)
-
-myNorm = function(in.dt,
-                  in.meas.col,
-                  in.rt.col = 'RealTime',
-                  in.rt.min = 10,
-                  in.rt.max = 20,
-                  in.by.cols = NULL,
-                  in.robust = TRUE,
-                  in.type = 'z.score') {
-  loc.dt <-
-    copy(in.dt) # copy so as not to alter original dt object w intermediate assignments
-  
-  if (is.null(in.by.cols)) {
-    if (in.robust)
-      loc.dt.pre.aggr = loc.dt[get(in.rt.col) >= in.rt.min &
-                                 get(in.rt.col) <= in.rt.max, .(meas.md = median(get(in.meas.col), na.rm = TRUE),
-                                                               meas.mad = mad(get(in.meas.col), na.rm = TRUE))]
-    else
-      loc.dt.pre.aggr = loc.dt[get(in.rt.col) >= in.rt.min &
-                                 get(in.rt.col) <= in.rt.max, .(meas.md = mean(get(in.meas.col), na.rm = TRUE),
-                                                               meas.mad = sd(get(in.meas.col), na.rm = TRUE))]
-    
-    loc.dt = cbind(loc.dt, loc.dt.pre.aggr)
-  }  else {
-    if (in.robust)
-      loc.dt.pre.aggr = loc.dt[get(in.rt.col) >= in.rt.min &
-                                 get(in.rt.col) <= in.rt.max, .(meas.md = median(get(in.meas.col), na.rm = TRUE),
-                                                               meas.mad = mad(get(in.meas.col), na.rm = TRUE)), by = in.by.cols]
-    else
-      loc.dt.pre.aggr = loc.dt[get(in.rt.col) >= in.rt.min &
-                                 get(in.rt.col) <= in.rt.max, .(meas.md = mean(get(in.meas.col), na.rm = TRUE),
-                                                               meas.mad = sd(get(in.meas.col), na.rm = TRUE)), by = in.by.cols]
-    
-    loc.dt = merge(loc.dt, loc.dt.pre.aggr, by = in.by.cols)
-  }
-  
-  
-  if (in.type == 'z.score') {
-    loc.dt[, meas.norm := (get(in.meas.col) - meas.md) / meas.mad]
-  } else {
-    loc.dt[, meas.norm := (get(in.meas.col) / meas.md)]
-  }
-  
-  setnames(loc.dt, 'meas.norm', paste0(in.meas.col, '.norm'))
-  
-  loc.dt[, c('meas.md', 'meas.mad') := NULL]
-  return(loc.dt)
-}
-
 # Plots a scatter plot with marginal histograms
 # Points are connected by a line (grouping by cellID)
 #

modules/tabAUC.R

@@ -6,10 +6,15 @@ modAUCplotUI =  function(id, label = "Plot Area Under Curves") {
   ns <- NS(id)
   
   tagList(
-      uiOutput(ns('uiSlTimeTrim')),
-      modStatsUI(ns('dispStats')),
-      br(),
-      modBoxPlotUI(ns('boxPlot')
+    h4(
+      "Calculate area under curve and plot per group"
+    ),
+    br(),
+    
+    uiOutput(ns('uiSlTimeTrim')),
+    modStatsUI(ns('dispStats')),
+    br(),
+    modBoxPlotUI(ns('boxPlot')
     )
   )
 }
@@ -35,23 +40,23 @@ modAUCplot = function(input, output, session, in.data, in.fname = 'boxplotAUC.pd
   output$uiSlTimeTrim = renderUI({
     cat(file = stderr(), 'UI uiSlTimeTrim\n')
     
-      locTpts  = getDataTpts()
-      
-      if(is.null(locTpts))
-        return(NULL)
-      
-      locRTmin = min(locTpts)
-      locRTmax = max(locTpts)
-      
-      sliderInput(
-        ns('slTimeTrim'),
-        label = 'Select time range for AUC calculation',
-        min = locRTmin,
-        max = locRTmax,
-        value = c(locRTmin, locRTmax),
-        step = 1
-      )
-      
+    locTpts  = getDataTpts()
+    
+    if(is.null(locTpts))
+      return(NULL)
+    
+    locRTmin = min(locTpts)
+    locRTmax = max(locTpts)
+    
+    sliderInput(
+      ns('slTimeTrim'),
+      label = 'Select time range for AUC calculation',
+      min = locRTmin,
+      max = locRTmax,
+      value = c(locRTmin, locRTmax),
+      step = 1
+    )
+    
   })
   
   AUCcells = reactive({
@@ -65,7 +70,7 @@ modAUCplot = function(input, output, session, in.data, in.fname = 'boxplotAUC.pd
       return(loc.res)
     }
   })
-
+  
   callModule(modStats, 'dispStats',
              in.data = AUCcells,
              in.meascol = 'AUC',

modules/tabClHier.R

@@ -148,12 +148,12 @@ clustHierUI <- function(id, label = "Hierarchical CLustering") {
                plotOutput(ns('outPlotHier'))
       ),
 
-      tabPanel('Time-courses',
-               modTrajPlotUI(ns('modPlotHierTraj'))),
-      
       tabPanel('Averages',
                modTrajRibbonPlotUI(ns('modPlotHierTrajRibbon'))),
       
+      tabPanel('Time-courses',
+               modTrajPlotUI(ns('modPlotHierTraj'))),
+      
       tabPanel('Cluster dist.',
                modClDistPlotUI(ns('hierClDistPlot'), 'xxx'))
       

modules/tabClHierSpar.R

@@ -159,8 +159,12 @@ clustHierSparUI <- function(id, label = "Sparse Hierarchical CLustering") {
                plotOutput(ns('outPlotHierSpar'))
       ),
 
+      tabPanel('Averages',
+               modTrajRibbonPlotUI(ns('modPlotHierSparTrajRibbon'))),
+      
       tabPanel('Time-courses',
                modTrajPlotUI(ns('modPlotHierSparTraj'))),
+      
       tabPanel('Cluster dist.',
                modClDistPlotUI(ns('hierClSparDistPlot')))
     )
@@ -222,9 +226,9 @@ clustHierSpar <- function(input, output, session, in.data4clust, in.data4trajPlo
       return()
     }
     
+    #cat('rownames: ', rownames(dm.t), '\n')
+    #cat('=============\ndimensions:', dim(dm.t), '\n')
     
-    cat('rownames: ', rownames(dm.t), '\n')
-    cat('=============\ndimensions:', dim(dm.t), '\n')
     perm.out <- HierarchicalSparseCluster.permute(
       dm.t,
       wbounds = NULL,
@@ -240,8 +244,8 @@ clustHierSpar <- function(input, output, session, in.data4clust, in.data4trajPlo
       dissimilarity = s.cl.spar.diss[as.numeric(input$selectPlotHierSparDiss)]
     )
     
-    cat('=============\nsparsehc:\n')
-    print(sparsehc$hc)
+    #cat('=============\nsparsehc:\n')
+    #print(sparsehc$hc)
     
     return(sparsehc)
   })
@@ -257,7 +261,7 @@ clustHierSpar <- function(input, output, session, in.data4clust, in.data4trajPlo
     
     dend <- as.dendrogram(sparsehc$hc)
     
-    cat('=============\ncutree:\n', dendextend::cutree(dend, input$inPlotHierSparNclust, order_clusters_as_data = TRUE), '\n')
+    #cat('=============\ncutree:\n', dendextend::cutree(dend, input$inPlotHierSparNclust, order_clusters_as_data = TRUE), '\n')
     
     dend <- color_branches(dend, 
                            col = rainbow_hcl,
@@ -322,7 +326,7 @@ clustHierSpar <- function(input, output, session, in.data4clust, in.data4trajPlo
     
     cat(file = stderr(), 'data4trajPlotClSpar: dt not NULL\n')
     
-    cat('rownames: ', rownames(in.data4clust()), '\n')
+    #cat('rownames: ', rownames(in.data4clust()), '\n')
 
     # get cellIDs with cluster assignments based on dendrogram cut
     loc.dt.cl = getDataClSpar(userFitDendHierSpar(), input$inPlotHierSparNclust, getDataTrackObjLabUni_afterTrim())

modules/trajRibbonPlot.R

@@ -45,7 +45,11 @@ modTrajRibbonPlotUI =  function(id, label = "Plot Individual Time Series") {
 }
 
 
-modTrajRibbonPlot = function(input, output, session, in.data, in.facet = 'group', in.facet.color = NULL, in.fname = 'tCourses.pdf') {
+modTrajRibbonPlot = function(input, output, session, 
+                             in.data, 
+                             in.facet = 'group', 
+                             in.facet.color = NULL, 
+                             in.fname = 'tCoursesMeans.pdf') {
   
   ns <- session$ns
   
@@ -166,7 +170,7 @@ modTrajRibbonPlot = function(input, output, session, in.data, in.facet = 'group'
       loc.facet.col = loc.facet.col[loc.groups]
     }
     
-    loc.dt.aggr = calcTrajCI(in.dt = loc.dt, 
+    loc.dt.aggr = tca::calcTrajCI(in.dt = loc.dt, 
                              in.col.meas = 'y', 
                              in.col.by = c(in.facet, 'realtime'), 
                              in.type = 'normal')
@@ -181,22 +185,6 @@ modTrajRibbonPlot = function(input, output, session, in.data, in.facet = 'group'
                            xlab.arg = 'Time (min)',
                            ylab.arg = '')
     
-    # p.out = myGgplotTraj(
-    #   dt.arg = loc.dt,
-    #   x.arg = 'realtime',
-    #   y.arg = 'y',
-    #   group.arg = "id",
-    #   facet.arg = in.facet,
-    #   facet.ncol.arg = input$inPlotTrajFacetNcol,
-    #   facet.color.arg = loc.facet.col,
-    #   xlab.arg = 'Time (min)',
-    #   line.col.arg = loc.line.col.arg,
-    #   aux.label1 = if (locPos) 'pos.x' else NULL,
-    #   aux.label2 = if (locPos) 'pos.y' else NULL,
-    #   aux.label3 = if (locObjNum) 'obj.num' else NULL,
-    #   stat.arg = input$chBPlotTrajStat
-    # )
-    
     return(p.out)
   }
 }
\ No newline at end of file

server.R

@@ -21,6 +21,7 @@ library(sparcl) # sparse hierarchical and k-means
 library(scales) # for percentages on y scale
 library(dtw) # for dynamic time warping
 library(imputeTS) # for interpolating NAs
+library(tca) # for time series manipulatiom, e.g. normTraj, genTraj, plotTrajRibbon
 
 # increase file upload limit
 options(shiny.maxRequestSize = 200 * 1024 ^ 2)
@@ -53,7 +54,7 @@ shinyServer(function(input, output, session) {
   dataGen1 <- eventReactive(input$inDataGen1, {
     cat("dataGen1\n")
     
-    return(userDataGen())
+    return(tca::genTraj(in.nwells = 3))
   })
   
   # load main data file
@@ -171,7 +172,7 @@ shinyServer(function(input, output, session) {
   output$varSelSite = renderUI({
     cat(file = stderr(), 'UI varSelSite\n')
     
-    if (!input$chBtrackUni) {
+    if (input$chBtrackUni) {
       locCols = getDataNucCols()
       locColSel = locCols[grep('(S|s)ite|(S|s)eries', locCols)[1]] # index 1 at the end in case more matches; select 1st
       
@@ -339,13 +340,6 @@ shinyServer(function(input, output, session) {
   ####
   ## data processing
   
-  # generate random dataset 1
-  dataGen1 <- eventReactive(input$inDataGen1, {
-    cat("dataGen1\n")
-    
-    return(userDataGen())
-  })
-  
   dataInBoth <- reactive({
     # Without direct references to inDataGen1,2 and inFileLoad, inDataGen2
     #    does not trigger running this reactive once inDataGen1 is used.
@@ -412,7 +406,7 @@ shinyServer(function(input, output, session) {
     if (is.null(loc.dt))
       return(NULL)
     
-    if (!input$chBtrackUni) {
+    if (input$chBtrackUni) {
       loc.types = lapply(loc.dt, class)
       if(loc.types[[input$inSelTrackLabel]] %in% c('numeric', 'integer') & loc.types[[input$inSelSite]] %in% c('numeric', 'integer'))
       {
@@ -645,7 +639,7 @@ shinyServer(function(input, output, session) {
     ## Normalization
     # F-n myNorm adds additional column with .norm suffix
     if (input$chBnorm) {
-      loc.out = myNorm(
+      loc.out = tca::normTraj(
         in.dt = loc.out,
         in.meas.col = 'y',
         in.rt.col = 'realtime',
@@ -730,9 +724,16 @@ shinyServer(function(input, output, session) {
     }
   )
   
+  ###### Trajectory plotting
+  callModule(modTrajRibbonPlot, 'modTrajRibbon', data4trajPlot)
+  callModule(modTrajRibbonPlot, 'modTrajRibbon', 
+             in.data = data4trajPlot)
   
-  ####
-  ## UI for trajectory plot
+  ###### Trajectory plotting
+  callModule(modTrajPlot, 'modTrajPlot', data4trajPlot)
+  
+  ## UI for selecting trajectories
+  # The output data table of data4trajPlot is modified based on inSelHighlight field
   output$varSelHighlight = renderUI({
     cat(file = stderr(), 'UI varSelHighlight\n')
     
@@ -752,17 +753,12 @@ shinyServer(function(input, output, session) {
     }
   })
   
-  ###### Trajectory plotting
-  callModule(modTrajPlot, 'modTrajPlot', data4trajPlot)
-  
-  ###### AUC caluclation and plotting
+  ###### AUC calculation and plotting
   callModule(modAUCplot, 'tabAUC', data4trajPlot)
   
   ###### Box-plot
   callModule(tabBoxPlot, 'tabBoxPlot', data4trajPlot)
   
-  
-  
   ###### Scatter plot
   callModule(tabScatterPlot, 'tabScatter', data4trajPlot)
   

ui.R

@@ -42,7 +42,7 @@ shinyUI(fluidPage(
       uiOutput('uiButLoadTrajRem'),
       
       tags$hr(),
-      checkboxInput('chBtrackUni', 'Track Label unique across entire dataset', TRUE),
+      checkboxInput('chBtrackUni', 'Create unique TrackLabel', T),
       helpPopup(
         title = 'Create unique cell ID',
         content = help.text[2],
@@ -51,7 +51,7 @@ shinyUI(fluidPage(
       ),
       uiOutput('varSelSite'),
       uiOutput('varSelTrackLabel'),
-
+      
       tags$hr(),
       checkboxInput('chBgroup', 'Dataset contains grouping column (e.g. treatment, condition)', TRUE),                
       uiOutput('varSelGroup'),
@@ -87,47 +87,58 @@ shinyUI(fluidPage(
       downloadButton('downloadDataClean', 'Download mod\'d data')
     ),
     
-    mainPanel(tabsetPanel(
-      tabPanel(
-        "Time courses",
-        h4(
-          "Plot time series"
+    mainPanel(
+      tabsetPanel(
+        tabPanel(
+          "Time series",
+          h4(
+            "Plot time series: means per group or individual"
+          ),
+          br(),
+          
+          tabsetPanel(
+            tabPanel("Means",
+                     br(),
+                     modTrajRibbonPlotUI('modTrajRibbon')
+            ),
+            
+            tabPanel(
+              "Individual",
+              br(),
+              checkboxInput('chBhighlightTraj', 'Highlight trajectories?', FALSE),
+              uiOutput('varSelHighlight'),
+              br(),
+              modTrajPlotUI('modTrajPlot')
+            )
+          )
         ),
         
-        br(),
-        checkboxInput('chBhighlightTraj', 'Highlight trajectories?', FALSE),
-        uiOutput('varSelHighlight'),
-        br(),
-        modTrajPlotUI('modTrajPlot')
-      ),
-      
-      
-      tabPanel(
-        "AUC",
-        modAUCplotUI('tabAUC')
-      ),
-      
-      tabPanel(
-        "Box-plots",
-        tabBoxPlotUI('tabBoxPlot')
-      ),
-      
-      
-      # scatter plot
-      tabPanel(
-        'Scatter',
-        tabScatterPlotUI('tabScatter')
-      ),
-      
-      tabPanel(
-        'Hierarchical',
-        clustHierUI('tabClHier')
-      ),
-      
-      tabPanel(
-        'Hierarchical Sparse',
-        clustHierSparUI('tabClHierSpar')
-      )
-    ))
+        tabPanel(
+          "AUC",
+          modAUCplotUI('tabAUC')
+        ),
+        
+        tabPanel(
+          "Box-plots",
+          tabBoxPlotUI('tabBoxPlot')
+        ),
+        
+        
+        # scatter plot
+        tabPanel(
+          'Scatter',
+          tabScatterPlotUI('tabScatter')
+        ),
+        
+        tabPanel(
+          'Hierarchical',
+          clustHierUI('tabClHier')
+        ),
+        
+        tabPanel(
+          'Hierarchical Sparse',
+          clustHierSparUI('tabClHierSpar')
+        )
+      ))
   )
 ))
\ No newline at end of file