server.R 26.5 KB
Newer Older
dmattek's avatar
dmattek committed
1
#
dmattek's avatar
dmattek committed
2
3
4
5
# Time Course Inspector: Shiny app for plotting time series data
# Author: Maciej Dobrzynski
#
# This is the server logic for a Shiny web application.
dmattek's avatar
dmattek committed
6
7
8
9
10
11
#

library(shiny)
library(shinyjs) #http://deanattali.com/shinyjs/
library(data.table)
library(ggplot2)
dmattek's avatar
dmattek committed
12
library(gplots) # for heatmap.2
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
13
14
15
library(plotly) # interactive plot
library(DT) # interactive tables

dmattek's avatar
dmattek committed
16
library(dendextend) # for color_branches
17
library(colorspace) # for palettes (used to colour dendrogram)
dmattek's avatar
dmattek committed
18
19
library(RColorBrewer)
library(scales) # for percentages on y scale
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
20
21

library(sparcl) # sparse hierarchical and k-means
dmattek's avatar
Added:    
dmattek committed
22
23
library(dtw) # for dynamic time warping
library(imputeTS) # for interpolating NAs
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
24
25
26
27
library(robust) # for robust linear regression
library(MASS)
library(pracma) # for trapz

dmattek's avatar
dmattek committed
28

29
# Global parameters ----
dmattek's avatar
dmattek committed
30
# change to increase the limit of the upload file size
dmattek's avatar
Added:    
dmattek committed
31
options(shiny.maxRequestSize = 200 * 1024 ^ 2)
dmattek's avatar
dmattek committed
32

dmattek's avatar
dmattek committed
33
# Server logic ----
dmattek's avatar
dmattek committed
34
shinyServer(function(input, output, session) {
35
  useShinyjs()
dmattek's avatar
dmattek committed
36
  
37
  # This is only set at session start
dmattek's avatar
dmattek committed
38
  # We use this as a way to determine which input was
39
40
  # clicked in the dataInBoth reactive
  counter <- reactiveValues(
dmattek's avatar
dmattek committed
41
42
43
    # The value of actionButton is the number of times the button is pressed
    dataGen1        = isolate(input$inDataGen1),
    dataLoadNuc     = isolate(input$inButLoadNuc),
44
45
    dataLoadTrajRem = isolate(input$inButLoadTrajRem),
    dataLoadStim    = isolate(input$inButLoadStim)
dmattek's avatar
dmattek committed
46
  )
dmattek's avatar
dmattek committed
47
48
49
50
51
52
53
54
55

  nCellsCounter <- reactiveValues(
    nCellsOrig = 0,
    nCellsAfterOutlierTrim = 0
  )
    
  myReactVals = reactiveValues(
    outlierIDs = NULL
  )
dmattek's avatar
dmattek committed
56
  
dmattek's avatar
dmattek committed
57
  # UI-side-panel-data-load ----
dmattek's avatar
dmattek committed
58
  
dmattek's avatar
dmattek committed
59
  # Generate random dataset
60
  dataGen1 <- eventReactive(input$inDataGen1, {
61
    if (DEB)
62
      cat("server:dataGen1\n")
63
    
dmattek's avatar
dmattek committed
64
    return(LOCgenTraj(in.nwells = 3, in.addout = 3))
65
66
  })
  
dmattek's avatar
dmattek committed
67
  # Load main data file
68
  dataLoadNuc <- eventReactive(input$inButLoadNuc, {
69
    if (DEB)
70
      cat("server:dataLoadNuc\n")
71

72
73
74
75
76
77
78
    locFilePath = input$inFileLoadNuc$datapath
    
    counter$dataLoadNuc <- input$inButLoadNuc - 1
    
    if (is.null(locFilePath) || locFilePath == '')
      return(NULL)
    else {
79
      return(fread(locFilePath, strip.white = T))
80
81
82
    }
  })
  
dmattek's avatar
dmattek committed
83
84
85
86
  # This button will reset the inFileLoad
  observeEvent(input$butReset, {
    reset("inFileLoadNuc")  # reset is a shinyjs function
  })
87

dmattek's avatar
dmattek committed
88
  # Load data with trajectories to remove
89
  dataLoadTrajRem <- eventReactive(input$inButLoadTrajRem, {
90
    if (DEB)
91
      cat(file = stdout(), "server:dataLoadTrajRem\n")
92
    
93
94
95
96
97
98
99
100
101
102
    locFilePath = input$inFileLoadTrajRem$datapath
    
    counter$dataLoadTrajRem <- input$inButLoadTrajRem - 1
    
    if (is.null(locFilePath) || locFilePath == '')
      return(NULL)
    else {
      return(fread(locFilePath))
    }
  })
dmattek's avatar
dmattek committed
103
  
dmattek's avatar
dmattek committed
104
  # Load data with stimulation pattern
105
  dataLoadStim <- eventReactive(input$inButLoadStim, {
106
    if (DEB)
107
      cat(file = stdout(), "server:dataLoadStim\n")
108
    
109
110
111
112
113
114
115
116
117
118
119
120
    locFilePath = input$inFileLoadStim$datapath
    
    counter$dataLoadStim <- input$inButLoadStim - 1
    
    if (is.null(locFilePath) || locFilePath == '')
      return(NULL)
    else {
      return(fread(locFilePath))
    }
  })
  
    
dmattek's avatar
Added:    
dmattek committed
121
122
  # UI for loading csv with cell IDs for trajectory removal
  output$uiFileLoadTrajRem = renderUI({
123
    if (DEB)
124
      cat(file = stdout(), 'server:uiFileLoadTrajRem\n')
dmattek's avatar
Added:    
dmattek committed
125
126
127
128
    
    if(input$chBtrajRem) 
      fileInput(
        'inFileLoadTrajRem',
129
        'Select file and press "Load Data"',
dmattek's avatar
Added:    
dmattek committed
130
131
132
133
134
        accept = c('text/csv', 'text/comma-separated-values,text/plain')
      )
  })
  
  output$uiButLoadTrajRem = renderUI({
135
    if (DEB)
136
      cat(file = stdout(), 'server:uiButLoadTrajRem\n')
dmattek's avatar
Added:    
dmattek committed
137
138
139
140
141
    
    if(input$chBtrajRem)
      actionButton("inButLoadTrajRem", "Load Data")
  })

142
143
  # UI for loading csv with stimulation pattern
  output$uiFileLoadStim = renderUI({
144
    if (DEB)
145
      cat(file = stdout(), 'server:uiFileLoadStim\n')
dmattek's avatar
Added:    
dmattek committed
146
    
147
148
149
    if(input$chBstim) 
      fileInput(
        'inFileLoadStim',
150
        'Select file and press "Load Data"',
151
152
153
154
155
        accept = c('text/csv', 'text/comma-separated-values,text/plain')
      )
  })
  
  output$uiButLoadStim = renderUI({
156
    if (DEB)
157
      cat(file = stdout(), 'server:uiButLoadStim\n')
dmattek's avatar
Added:    
dmattek committed
158
    
159
160
    if(input$chBstim)
      actionButton("inButLoadStim", "Load Data")
dmattek's avatar
Added:    
dmattek committed
161
162
  })
  
163

dmattek's avatar
dmattek committed
164
  
dmattek's avatar
dmattek committed
165
  # UI-side-panel-column-selection ----
dmattek's avatar
dmattek committed
166
  output$varSelTrackLabel = renderUI({
167
    if (DEB)
168
      cat(file = stdout(), 'server:varSelTrackLabel\n')
169
    
dmattek's avatar
dmattek committed
170
    locCols = getDataNucCols()
171
    locColSel = locCols[grep('(T|t)rack|ID|id', locCols)[1]] # index 1 at the end in case more matches; select 1st; matches TrackLabel, tracklabel, Track Label etc
dmattek's avatar
dmattek committed
172
173
174
    
    selectInput(
      'inSelTrackLabel',
dmattek's avatar
dmattek committed
175
      'Select ID:',
dmattek's avatar
dmattek committed
176
177
178
179
180
181
182
      locCols,
      width = '100%',
      selected = locColSel
    )
  })
  
  output$varSelTime = renderUI({
183
    if (DEB)
184
      cat(file = stdout(), 'server:varSelTime\n')
185
    
dmattek's avatar
dmattek committed
186
    locCols = getDataNucCols()
187
    locColSel = locCols[grep('(T|t)ime|Metadata_T', locCols)[1]] # index 1 at the end in case more matches; select 1st; matches RealTime, realtime, real time, etc.
dmattek's avatar
dmattek committed
188
189
190
    
    selectInput(
      'inSelTime',
191
      'Select time column:',
dmattek's avatar
dmattek committed
192
193
194
195
196
      locCols,
      width = '100%',
      selected = locColSel
    )
  })
197
198

  output$varSelTimeFreq = renderUI({
199
    if (DEB)
200
      cat(file = stdout(), 'server:varSelTimeFreq\n')
201
    
202
203
204
    if (input$chBtrajInter) {
      numericInput(
        'inSelTimeFreq',
205
        'Frequency of time units:',
206
207
208
209
210
211
        min = 1,
        step = 1,
        width = '100%',
        value = 1
      )
    }
212
  })
dmattek's avatar
dmattek committed
213
  
dmattek's avatar
dmattek committed
214
  # This is the main field to select plot facet grouping
dmattek's avatar
dmattek committed
215
  # It's typically a column with the entire experimental description,
dmattek's avatar
dmattek committed
216
217
  # e.g.1 Stim_All_Ch or Stim_All_S.
  # e.g.2 a combination of 3 columns called Stimulation_...
dmattek's avatar
dmattek committed
218
  output$varSelGroup = renderUI({
219
    if (DEB)
220
      cat(file = stdout(), 'server:varSelGroup\n')
dmattek's avatar
dmattek committed
221
    
dmattek's avatar
dmattek committed
222
223
224
225
226
    if (input$chBgroup) {
      
      locCols = getDataNucCols()
      
      if (!is.null(locCols)) {
227
        locColSel = locCols[grep('(G|g)roup|(S|s)tim|(S|s)timulation|(S|s)ite', locCols)[1]]
228
229

        #cat('UI varSelGroup::locColSel ', locColSel, '\n')
dmattek's avatar
dmattek committed
230
231
        selectInput(
          'inSelGroup',
232
          'Select columns for plot grouping:',
dmattek's avatar
dmattek committed
233
234
235
236
237
          locCols,
          width = '100%',
          selected = locColSel,
          multiple = TRUE
        )
dmattek's avatar
dmattek committed
238
239
240
241
      }
    }
  })
  
242
243
  # UI for selecting grouping to add to track ID to make 
  # the track ID unique across entire dataset
dmattek's avatar
dmattek committed
244
  output$varSelSite = renderUI({
245
    if (DEB)
246
      cat(file = stdout(), 'server:varSelSite\n')
dmattek's avatar
dmattek committed
247
    
248
    if (input$chBtrackUni) {
dmattek's avatar
Added:    
dmattek committed
249
      locCols = getDataNucCols()
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
250
      locColSel = locCols[grep('(S|s)ite|(S|s)eries|(F|f)ov|(G|g)roup', locCols)[1]] # index 1 at the end in case more matches; select 1st
dmattek's avatar
Added:    
dmattek committed
251
252
253
      
      selectInput(
        'inSelSite',
254
        'Select grouping columns to add to track label:',
dmattek's avatar
Added:    
dmattek committed
255
256
        locCols,
        width = '100%',
257
258
        selected = locColSel,
        multiple = T
dmattek's avatar
Added:    
dmattek committed
259
260
      )
    }
dmattek's avatar
dmattek committed
261
262
263
264
  })
  
  
  output$varSelMeas1 = renderUI({
265
    if (DEB)
266
      cat(file = stdout(), 'server:varSelMeas1\n')
dmattek's avatar
dmattek committed
267
268
269
    locCols = getDataNucCols()
    
    if (!is.null(locCols)) {
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
270
      locColSel = locCols[grep('(R|r)atio|(I|i)ntensity|(Y|y)|(M|m)eas', locCols)[1]] # index 1 at the end in case more matches; select 1st
dmattek's avatar
dmattek committed
271

dmattek's avatar
dmattek committed
272
273
      selectInput(
        'inSelMeas1',
274
        'Select 1st meas.:',
dmattek's avatar
dmattek committed
275
276
277
278
279
280
281
282
283
        locCols,
        width = '100%',
        selected = locColSel
      )
    }
  })
  
  
  output$varSelMeas2 = renderUI({
284
    if (DEB)
285
      cat(file = stdout(), 'server:varSelMeas2\n')
286
    
dmattek's avatar
dmattek committed
287
288
289
290
    locCols = getDataNucCols()
    
    if (!is.null(locCols) &&
        !(input$inSelMath %in% c('', '1 / '))) {
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
291
      locColSel = locCols[grep('(R|r)atio|(I|i)ntensity|(Y|y)|(M|m)eas', locCols)[1]] # index 1 at the end in case more matches; select 1st
dmattek's avatar
dmattek committed
292

dmattek's avatar
dmattek committed
293
294
      selectInput(
        'inSelMeas2',
295
        'Select 2nd meas.',
dmattek's avatar
dmattek committed
296
297
298
299
300
301
302
        locCols,
        width = '100%',
        selected = locColSel
      )
    }
  })
  
dmattek's avatar
dmattek committed
303
  # UI-side-panel-trim x-axis (time) ----
dmattek's avatar
dmattek committed
304
  output$uiSlTimeTrim = renderUI({
305
    if (DEB)
306
      cat(file = stdout(), 'server:uiSlTimeTrim\n')
dmattek's avatar
dmattek committed
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
    
    if (input$chBtimeTrim) {
      locTpts  = getDataTpts()
      
      if(is.null(locTpts))
        return(NULL)
      
      locRTmin = min(locTpts)
      locRTmax = max(locTpts)
      
      sliderInput(
        'slTimeTrim',
        label = 'Time range to include',
        min = locRTmin,
        max = locRTmax,
        value = c(locRTmin, locRTmax),
        step = 1
      )
      
    }
  })
dmattek's avatar
dmattek committed
328
  
dmattek's avatar
dmattek committed
329
  # UI-side-panel-normalization ----
330
331
332
333
  
  # select normalisation method
  # - fold-change calculates fold change with respect to the mean
  # - z-score calculates z-score of the selected regione of the time series
dmattek's avatar
dmattek committed
334
  output$uiChBnorm = renderUI({
335
    if (DEB)
336
      cat(file = stdout(), 'server:uiChBnorm\n')
dmattek's avatar
dmattek committed
337
338
    
    if (input$chBnorm) {
339
      tagList(
dmattek's avatar
dmattek committed
340
341
342
      radioButtons(
        'rBnormMeth',
        label = 'Select method',
343
344
345
346
        choices = list('fold-change' = 'mean', 'z-score' = 'z.score'),
        width = "40%"
      ),
      bsTooltip('rBnormMeth', help.text.short[11], placement = "right", trigger = "hover", options = NULL)
dmattek's avatar
dmattek committed
347
348
349
350
      )
    }
  })
  
351
  # select the region of the time series for normalisation
dmattek's avatar
dmattek committed
352
  output$uiSlNorm = renderUI({
353
    if (DEB)
354
      cat(file = stdout(), 'server:uiSlNorm\n')
dmattek's avatar
dmattek committed
355
356
357
358
359
360
361
362
363
364
    
    if (input$chBnorm) {
      locTpts  = getDataTpts()
      
      if(is.null(locTpts))
        return(NULL)
      
      locRTmin = min(locTpts)
      locRTmax = max(locTpts)
      
365
      tagList(
dmattek's avatar
dmattek committed
366
367
      sliderInput(
        'slNormRtMinMax',
368
        label = 'Time span',
dmattek's avatar
dmattek committed
369
370
        min = locRTmin,
        max = locRTmax,
dmattek's avatar
dmattek committed
371
372
        value = c(locRTmin, 0.1 * locRTmax), 
        step = 1
373
374
      ),
      bsTooltip('slNormRtMinMax', help.text.short[12], placement = "right", trigger = "hover", options = NULL)
dmattek's avatar
dmattek committed
375
376
377
378
      )
    }
  })
  
379
  # use robust stats (median instead of mean, mad instead of sd)
dmattek's avatar
dmattek committed
380
  output$uiChBnormRobust = renderUI({
381
    if (DEB)
382
      cat(file = stdout(), 'server:uiChBnormRobust\n')
dmattek's avatar
dmattek committed
383
384
    
    if (input$chBnorm) {
385
      tagList(
dmattek's avatar
dmattek committed
386
387
      checkboxInput('chBnormRobust',
                    label = 'Robust stats',
388
389
390
391
                    FALSE, 
                    width = "40%"),
      bsTooltip('chBnormRobust', help.text.short[13], placement = "right", trigger = "hover", options = NULL)
      )
dmattek's avatar
dmattek committed
392
393
394
    }
  })
  
395
  # choose whether normalisation should be calculated for the entire dataset, group, or trajectory
dmattek's avatar
dmattek committed
396
  output$uiChBnormGroup = renderUI({
397
    if (DEB)
398
      cat(file = stdout(), 'server:uiChBnormGroup\n')
dmattek's avatar
dmattek committed
399
400
    
    if (input$chBnorm) {
401
      tagList(
dmattek's avatar
dmattek committed
402
      radioButtons('chBnormGroup',
dmattek's avatar
Mod:    
dmattek committed
403
                   label = 'Normalisation grouping',
404
405
406
407
                   choices = list('Entire dataset' = 'none', 'Per group' = 'group', 'Per trajectory' = 'id'), 
                   width = "40%"),
      bsTooltip('chBnormGroup', help.text.short[14], placement = "right", trigger = "hover", options = NULL)
      )
dmattek's avatar
dmattek committed
408
409
410
411
    }
  })
  
  
dmattek's avatar
dmattek committed
412
  
dmattek's avatar
dmattek committed
413

dmattek's avatar
dmattek committed
414
  # Processing-data ----
dmattek's avatar
dmattek committed
415
  
416
417
418
419
420
421
422
423
424
425
426
427
  dataInBoth <- reactive({
    # Without direct references to inDataGen1,2 and inFileLoad, inDataGen2
    #    does not trigger running this reactive once inDataGen1 is used.
    # This is one of the more nuanced areas of reactive programming in shiny
    #    due to the if else logic, it isn't fetched once inDataGen1 is available
    # The morale is use direct retrieval of inputs to guarantee they are available
    #    for if else logic checks!
    
    locInGen1 = input$inDataGen1
    locInLoadNuc = input$inButLoadNuc
    #locInLoadStim = input$inButLoadStim
    
428
    # Don't wrap around if(DEB) !!!
429
    cat(
430
      "server:dataInBoth\n   inGen1: ",
431
      locInGen1,
432
      "      prev=",
433
      isolate(counter$dataGen1),
434
      "\n   inDataNuc: ",
435
436
437
438
439
440
441
442
443
444
      locInLoadNuc,
      "   prev=",
      isolate(counter$dataLoadNuc),
      # "\ninDataStim: ",
      # locInLoadStim,
      # "   prev=",
      # isolate(counter$dataLoadStim),
      "\n"
    )
    
445
    # isolate the checks of the counter reactiveValues
446
447
    # as we set the values in this same reactive
    if (locInGen1 != isolate(counter$dataGen1)) {
448
      cat("server:dataInBoth if inDataGen1\n")
449
450
451
452
      dm = dataGen1()
      # no need to isolate updating the counter reactive values!
      counter$dataGen1 <- locInGen1
    } else if (locInLoadNuc != isolate(counter$dataLoadNuc)) {
453
      cat("server:dataInBoth if inDataLoadNuc\n")
454
      dm = dataLoadNuc()
455
456
457
458
459
460
461
462
463
464
465
466
467
      
      # convert to long format if radio box set to "wide"
      # the input data in long format should contain:
      # - the first row with a header: ID, 1, 2, 3...
      # - consecutive rows with time series, where columns are time points
      if (input$inRbutLongWide == 1) {
        # long to wide
        dm = melt(dm, id.vars = names(dm)[1], variable.name = COLRT, value.name = COLY)

        # convert column names with time points to a number
        dm[, (COLRT) := as.numeric(levels(get(COLRT)))[get(COLRT)]]
      }
      
468
469
470
      # no need to isolate updating the counter reactive values!
      counter$dataLoadNuc <- locInLoadNuc
    } else {
471
      cat("server:dataInBoth else\n")
472
473
474
475
476
477
      dm = NULL
    }
    return(dm)
  })
  
  # return column names of the main dt
dmattek's avatar
dmattek committed
478
  getDataNucCols <- reactive({
479
    if (DEB)
480
      cat(file = stdout(), 'server:getDataNucCols: in\n')
481
    
482
483
484
485
486
487
488
489
490
491
    loc.dt = dataInBoth()
    
    if (is.null(loc.dt))
      return(NULL)
    else
      return(colnames(loc.dt))
  })
  
  # return dt with an added column with unique track object label
  dataMod <- reactive({
492
    if (DEB)
493
      cat(file = stdout(), 'server:dataMod\n')
494
    
495
496
    loc.dt = dataInBoth()
    
dmattek's avatar
dmattek committed
497
    if (is.null(loc.dt))
498
499
      return(NULL)
    
500
    if (input$chBtrackUni) {
501
      # create unique track ID based on columns specified in input$inSelSite field and combine with input$inSelTrackLabel
502
      loc.dt[, (COLIDUNI) := do.call(paste, c(.SD, sep = "_")), .SDcols = c(input$inSelSite, input$inSelTrackLabel) ]
dmattek's avatar
Added:    
dmattek committed
503
    } else {
504
      # stay with track ID provided in the loaded dataset; has to be unique
505
      loc.dt[, (COLIDUNI) := get(input$inSelTrackLabel)]
dmattek's avatar
Added:    
dmattek committed
506
507
    }
    
dmattek's avatar
dmattek committed
508
    
dmattek's avatar
Added:    
dmattek committed
509
510
511
    # remove trajectories based on uploaded csv

    if (input$chBtrajRem) {
512
      if (DEB)
513
        cat(file = stdout(), 'server:dataMod: trajRem not NULL\n')
dmattek's avatar
Added:    
dmattek committed
514
515
      
      loc.dt.rem = dataLoadTrajRem()
dmattek's avatar
dmattek committed
516
      loc.dt = loc.dt[!(trackObjectsLabelUni %in% loc.dt.rem[[1]])]
dmattek's avatar
Added:    
dmattek committed
517
518
    }
    
519
520
521
    return(loc.dt)
  })
  
dmattek's avatar
dmattek committed
522
523
524
  # return all unique track object labels (created in dataMod)
  # This will be used to display in UI for trajectory highlighting
  getDataTrackObjLabUni <- reactive({
525
    if (DEB)
526
      cat(file = stdout(), 'server:getDataTrackObjLabUni\n')
527
    
dmattek's avatar
dmattek committed
528
    loc.dt = dataMod()
529
    
dmattek's avatar
dmattek committed
530
531
532
533
    if (is.null(loc.dt))
      return(NULL)
    else
      return(unique(loc.dt$trackObjectsLabelUni))
534
535
  })
  
dmattek's avatar
Mod:    
dmattek committed
536
  
dmattek's avatar
dmattek committed
537
538
539
  # return all unique time points (real time)
  # This will be used to display in UI for box-plot
  # These timepoints are from the original dt and aren't affected by trimming of x-axis
dmattek's avatar
dmattek committed
540
  getDataTpts <- reactive({
541
    if (DEB)
542
      cat(file = stdout(), 'server:getDataTpts\n')
543
    
dmattek's avatar
dmattek committed
544
    loc.dt = dataMod()
545
    
dmattek's avatar
dmattek committed
546
547
548
549
    if (is.null(loc.dt))
      return(NULL)
    else
      return(unique(loc.dt[[input$inSelTime]]))
550
551
  })
  
dmattek's avatar
dmattek committed
552
  
553
554
555
  
  # prepare data for plotting time courses
  # returns dt with these columns:
dmattek's avatar
dmattek committed
556
  #    realtime - selected from input
dmattek's avatar
dmattek committed
557
  #    y        - measurement selected from input
dmattek's avatar
dmattek committed
558
  #               (can be a single column or result of an operation on two cols)
559
560
  #    id       - trackObjectsLabelUni; created in dataMod based on TrackObjects_Label
  #               and FOV column such as Series or Site (if TrackObjects_Label not unique across entire dataset)
dmattek's avatar
dmattek committed
561
562
  #    group    - grouping variable for facetting from input
  #    mid.in   - column with trajectory selection status from the input file or
563
564
565
566
  #               highlight status from UI 
  #               (column created if mid.in present in uploaded data or tracks are selected in the UI)
  #    obj.num  - created if ObjectNumber column present in the input data 
  #    pos.x,y  - created if columns with x and y positions present in the input data
567
  data4trajPlot <- reactive({
568
    if (DEB)
569
      cat(file = stdout(), 'server:data4trajPlot\n')
570
571
    
    loc.dt = dataMod()
dmattek's avatar
dmattek committed
572
    if (is.null(loc.dt))
573
574
      return(NULL)
    
575
    # create expression for 'y' column based on measurements and math operations selected in UI
dmattek's avatar
dmattek committed
576
    if (input$inSelMath == '')
577
578
579
580
581
582
      loc.s.y = input$inSelMeas1
    else if (input$inSelMath == '1 / ')
      loc.s.y = paste0(input$inSelMath, input$inSelMeas1)
    else
      loc.s.y = paste0(input$inSelMeas1, input$inSelMath, input$inSelMeas2)
    
583
    # create expression for 'group' column
584
585
    # creates a merged column based on other columns from input
    # used for grouping of plot facets
dmattek's avatar
dmattek committed
586
587
588
589
590
591
592
593
594
595
596
    if (input$chBgroup) {
      if(length(input$inSelGroup) == 0)
        return(NULL)
      
      loc.s.gr = sprintf("paste(%s, sep=';')",
                         paste(input$inSelGroup, sep = '', collapse = ','))
    } else {
      # if no grouping required, fill 'group' column with 0
      # because all the plotting relies on the presence of the group column
      loc.s.gr = "paste('0')"
    }
597
    
dmattek's avatar
dmattek committed
598
599

    # column name with time
600
601
    loc.s.rt = input$inSelTime
    
dmattek's avatar
dmattek committed
602
603
    # Assign tracks selected for highlighting in UI
    loc.tracks.highlight = input$inSelHighlight
604
    locButHighlight = input$chBhighlightTraj
dmattek's avatar
dmattek committed
605
    
dmattek's avatar
Added:    
dmattek committed
606
607
    
    # Find column names with position
608
    loc.s.pos.x = names(loc.dt)[grep('(L|l)ocation.*X|(P|p)os.x|(P|p)osx', names(loc.dt))[1]]
609
    loc.s.pos.y = names(loc.dt)[grep('(L|l)ocation.*Y|(P|p)os.y|(P|p)osy', names(loc.dt))[1]]
dmattek's avatar
Added:    
dmattek committed
610
    
611
    if (DEB)
612
      cat('server:data4trajPlot:\n   Position columns: ', loc.s.pos.x, loc.s.pos.y, '\n')
613
614
    
    if (!is.na(loc.s.pos.x) & !is.na(loc.s.pos.y))
dmattek's avatar
Added:    
dmattek committed
615
616
617
618
      locPos = TRUE
    else
      locPos = FALSE
    
619
620
621
622
    
    # Find column names with ObjectNumber
    # This is different from TrackObject_Label and is handy to keep
    # because labels on segmented images are typically ObjectNumber
623
624
625
626
627
628
    loc.s.objnum = names(loc.dt)[grep('(O|o)bject(N|n)umber', names(loc.dt))[1]]
    #cat('data4trajPlot::loc.s.objnum ', loc.s.objnum, '\n')
    if (is.na(loc.s.objnum)) {
      locObjNum = FALSE
    }
    else {
dmattek's avatar
dmattek committed
629
      loc.s.objnum = loc.s.objnum[1]
630
      locObjNum = TRUE
dmattek's avatar
dmattek committed
631
    }
632
633
    
    
634
635
    # if dataset contains column mid.in with trajectory filtering status,
    # then, include it in plotting
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
    if (sum(names(loc.dt) %in% 'mid.in') > 0)
      locMidIn = TRUE
    else
      locMidIn = FALSE
    
    ## Build expression for selecting columns from loc.dt
    # Core columns
    s.colexpr = paste0('.(y = ', loc.s.y,
                       ', id = trackObjectsLabelUni', 
                       ', group = ', loc.s.gr,
                       ', realtime = ', loc.s.rt)
    
    # account for the presence of 'mid.in' column in uploaded data
    if(locMidIn)
      s.colexpr = paste0(s.colexpr, 
                         ', mid.in = mid.in')
    
    # include position x,y columns in uploaded data
    if(locPos)
      s.colexpr = paste0(s.colexpr, 
                         ', pos.x = ', loc.s.pos.x,
                         ', pos.y = ', loc.s.pos.y)
    
    # include ObjectNumber column
    if(locObjNum)
      s.colexpr = paste0(s.colexpr, 
                         ', obj.num = ', loc.s.objnum)
    
    # close bracket, finish the expression
    s.colexpr = paste0(s.colexpr, ')')
    
    # create final dt for output based on columns selected above
    loc.out = loc.dt[, eval(parse(text = s.colexpr))]
    
670
671
672
673
674
    # Convert track ID to a factor.
    # This is necessary for, e.g. merging data with cluster assignments.
    # If input dataset has track ID as a number, such a merge would fail.
    loc.out[, (COLID) := as.factor(get(COLID))]
    
675
676
677
678
679
680
681
682
    
    # if track selection ON
    if (locButHighlight){
      # add a 3rd level with status of track selection
      # to a column with trajectory filtering status in the uploaded file
      if(locMidIn)
        loc.out[, mid.in := ifelse(id %in% loc.tracks.highlight, 'SELECTED', mid.in)]
      else
dmattek's avatar
Mod:    
dmattek committed
683
        # add a column with status of track selection
684
        loc.out[, mid.in := ifelse(id %in% loc.tracks.highlight, 'SELECTED', 'NOT SEL')]
685
    }
686
      
dmattek's avatar
dmattek committed
687

688
    ## Interpolate missing data and NA data points
689
    # From: https://stackoverflow.com/questions/28073752/r-how-to-add-rows-for-missing-values-for-unique-group-sequences
690
    # Tracks are interpolated only within first and last time points of every track id
691
692
    # Datasets can have different realtime frequency (e.g. every 1', 2', etc),
    # or the frame number metadata can be missing, as is the case for tCourseSelected files that already have realtime column.
693
    # Therefore, we cannot rely on that info to get time frequency; user must provide this number!
694
    
695
    setkeyv(loc.out, c(COLGR, COLID, COLRT))
696

697
698
    if (input$chBtrajInter) {
      # here we fill missing data with NA's
dmattek's avatar
dmattek committed
699
      loc.out = loc.out[setkeyv(loc.out[, .(seq(min(get(COLRT), na.rm = T), max(get(COLRT), na.rm = T), input$inSelTimeFreq)), by = c(COLGR, COLID)], c(COLGR, COLID, 'V1'))]
700
701
      
      # x-check: print all rows with NA's
702
      if (DEB) {
703
        cat(file = stdout(), 'server:data4trajPlot: Rows with NAs:\n')
704
705
        print(loc.out[rowSums(is.na(loc.out)) > 0, ])
      }
706
707
708
709
      
      # NA's may be already present in the dataset'.
      # Interpolate (linear) them with na.interpolate as well
      if(locPos)
dmattek's avatar
dmattek committed
710
        s.cols = c(COLY, COLPOSX, COLPOSY)
711
      else
dmattek's avatar
dmattek committed
712
        s.cols = c(COLY)
713
      
714
715
716
717
718
719
720
721
722
723
724
725
      # Interpolated columns should be of type numeric (float)
      # This is to ensure that interpolated columns are of porper type.
      
      # Apparently the loop is faster than lapply+SDcols
      for(col in s.cols) {
        #loc.out[, (col) := as.numeric(get(col))]
        data.table::set(loc.out, j = col, value = as.numeric(loc.out[[col]]))

        loc.out[, (col) := na.interpolation(get(col)), by = c(COLID)]        
      }
      
      # loc.out[, (s.cols) := lapply(.SD, na.interpolation), by = c(COLID), .SDcols = s.cols]
726
727
728
729
730
731
732
733
734
735
736
737
738
739
      
      
      # !!! Current issue with interpolation:
      # The column mid.in is not taken into account.
      # If a trajectory is selected in the UI,
      # the mid.in column is added (if it doesn't already exist in the dataset),
      # and for the interpolated point, it will still be NA. Not really an issue.
      #
      # Also, think about the current option of having mid.in column in the uploaded dataset.
      # Keep it? Expand it?
      # Create a UI filed for selecting the column with mid.in data.
      # What to do with that column during interpolation (see above)
      
    }    
dmattek's avatar
Mod:    
dmattek committed
740
    
741
    ## Trim x-axis (time)
dmattek's avatar
dmattek committed
742
    if(input$chBtimeTrim) {
dmattek's avatar
dmattek committed
743
      loc.out = loc.out[get(COLRT) >= input$slTimeTrim[[1]] & get(COLRT) <= input$slTimeTrim[[2]] ]
dmattek's avatar
dmattek committed
744
    }
dmattek's avatar
dmattek committed
745
    
746
    ## Normalization
dmattek's avatar
dmattek committed
747
    # F-n normTraj adds additional column with .norm suffix
dmattek's avatar
dmattek committed
748
    if (input$chBnorm) {
dmattek's avatar
dmattek committed
749
      loc.out = LOCnormTraj(
dmattek's avatar
dmattek committed
750
        in.dt = loc.out,
dmattek's avatar
dmattek committed
751
752
        in.meas.col = COLY,
        in.rt.col = COLRT,
dmattek's avatar
dmattek committed
753
754
755
756
757
758
759
        in.rt.min = input$slNormRtMinMax[1],
        in.rt.max = input$slNormRtMinMax[2],
        in.type = input$rBnormMeth,
        in.robust = input$chBnormRobust,
        in.by.cols = if(input$chBnormGroup %in% 'none') NULL else input$chBnormGroup
      )
      
dmattek's avatar
dmattek committed
760
761
      # Column with normalized data is renamed to the original name
      # Further code assumes column name y produced by data4trajPlot
Maciej Dobrzynski's avatar
Maciej Dobrzynski committed
762
763
      
      loc.out[, c(COLY) := NULL]
dmattek's avatar
dmattek committed
764
      setnames(loc.out, 'y.norm', COLY)
dmattek's avatar
dmattek committed
765
766
767
    }
    
    return(loc.out)
dmattek's avatar
dmattek committed
768
769
  })
  
dmattek's avatar
dmattek committed
770
771
772
773
774
775
  
  # prepare data for clustering
  # return a matrix with:
  # cells as columns
  # time points as rows
  data4clust <- reactive({
776
    if (DEB)  
777
      cat(file = stdout(), 'server:data4clust\n')
dmattek's avatar
dmattek committed
778
    
dmattek's avatar
dmattek committed
779
    loc.dt = data4trajPlotNoOut()
dmattek's avatar
dmattek committed
780
781
782
    if (is.null(loc.dt))
      return(NULL)
    
783
784
785
786
    # convert from long to wide format
    loc.dt.wide = dcast(loc.dt, 
                    reformulate(response = COLID, termlabels = COLRT), 
                    value.var = COLY)
dmattek's avatar
dmattek committed
787
    
788
789
    # store row names for later
    loc.rownames = loc.dt.wide[[COLID]]
dmattek's avatar
Mod:    
dmattek committed
790
    
791
792
    # omit first column that contains row names
    loc.m.out = as.matrix(loc.dt.wide[, -1])
dmattek's avatar
Added:    
dmattek committed
793
    
794
795
    # assign row names to the matrix
    rownames(loc.m.out) = loc.rownames
dmattek's avatar
Added:    
dmattek committed
796
    
797
    return(loc.m.out)
dmattek's avatar
Mod:    
dmattek committed
798
  }) 
dmattek's avatar
dmattek committed
799
  
dmattek's avatar
dmattek committed
800
  
801
802
803
  # prepare data with stimulation pattern
  # this dataset is displayed underneath of trajectory plot (modules/trajPlot.R) as geom_segment
  data4stimPlot <- reactive({
804
    if (DEB)  
805
      cat(file = stdout(), 'server:data4stimPlot\n')
806
807
    
    if (input$chBstim) {
808
      if (DEB)  
809
        cat(file = stdout(), 'server:data4stimPlot: stim not NULL\n')
810
811
812
813
      
      loc.dt.stim = dataLoadStim()
      return(loc.dt.stim)
    } else {
814
      if (DEB)  
815
        cat(file = stdout(), 'server:data4stimPlot: stim is NULL\n')
816
      
817
818
819
820
      return(NULL)
    }
  })
  
dmattek's avatar
Added:    
dmattek committed
821
822
823
  # download data as prepared for plotting
  # after all modification
  output$downloadDataClean <- downloadHandler(
dmattek's avatar
dmattek committed
824
    filename = FCSVTCCLEAN,
dmattek's avatar
Added:    
dmattek committed
825
    content = function(file) {
dmattek's avatar
dmattek committed
826
      write.csv(data4trajPlotNoOut(), file, row.names = FALSE)
dmattek's avatar
Added:    
dmattek committed
827
828
829
    }
  )
  
dmattek's avatar
dmattek committed
830
831
832
  # Plotting-trajectories ----

  # UI for selecting trajectories
833
  # The output data table of data4trajPlot is modified based on inSelHighlight field
dmattek's avatar
dmattek committed
834
  output$varSelHighlight = renderUI({
835
    if (DEB)  
836
      cat(file = stdout(), 'server:varSelHighlight\n')
dmattek's avatar
dmattek committed
837
    
dmattek's avatar
dmattek committed
838
839
840
    locBut = input$chBhighlightTraj
    if (!locBut)
      return(NULL)
dmattek's avatar
dmattek committed
841
    
dmattek's avatar
dmattek committed
842
    loc.v = getDataTrackObjLabUni()
dmattek's avatar
dmattek committed
843
    if (!is.null(loc.v)) {
844
      selectInput(
dmattek's avatar
dmattek committed
845
        'inSelHighlight',
846
        'Select one or more trajectories:',
dmattek's avatar
dmattek committed
847
        loc.v,
848
        width = '100%',
dmattek's avatar
dmattek committed
849
        multiple = TRUE
850
      )
dmattek's avatar
dmattek committed
851
852
853
    }
  })
  
dmattek's avatar
dmattek committed
854
855
856
  # Taking out outliers 
  data4trajPlotNoOut = callModule(modSelOutliers, 'returnOutlierIDs', data4trajPlot)
  
dmattek's avatar
dmattek committed
857
858
  # Trajectory plotting - ribbon
  callModule(modTrajRibbonPlot, 'modTrajRibbon', 
dmattek's avatar
dmattek committed
859
             in.data = data4trajPlotNoOut,
dmattek's avatar
dmattek committed
860
             in.data.stim = data4stimPlot,
dmattek's avatar
dmattek committed
861
             in.fname = function() return(FPDFTCMEAN))
dmattek's avatar
dmattek committed
862
  
dmattek's avatar
dmattek committed
863
  # Trajectory plotting - individual
dmattek's avatar
dmattek committed
864
  callModule(modTrajPlot, 'modTrajPlot', 
dmattek's avatar
dmattek committed
865
             in.data = data4trajPlotNoOut, 
dmattek's avatar
dmattek committed
866
             in.data.stim = data4stimPlot,
dmattek's avatar
dmattek committed
867
             in.fname = function() {return(FPDFTCSINGLE)})
dmattek's avatar
dmattek committed
868
  
869
870
871
872
873
  # Trajectory plotting - PSD
  callModule(modPSDPlot, 'modPSDPlot',
             in.data = data4trajPlotNoOut,
             in.fname = function() {return(FPDFTCPSD)})
  
dmattek's avatar
dmattek committed
874
875
  
  # Tabs ----
876
  ###### AUC calculation and plotting
dmattek's avatar
dmattek committed
877
  callModule(modAUCplot, 'tabAUC', data4trajPlotNoOut, in.fname = function() return(FPDFBOXAUC))
dmattek's avatar
Added:    
dmattek committed
878
  
dmattek's avatar
Added:    
dmattek committed
879
  ###### Box-plot
dmattek's avatar
dmattek committed
880
  callModule(tabBoxPlot, 'tabBoxPlot', data4trajPlotNoOut, in.fname = function() return(FPDFBOXTP))
dmattek's avatar
dmattek committed
881
  
dmattek's avatar
dmattek committed
882
  ###### Scatter plot
dmattek's avatar
dmattek committed
883
  callModule(tabScatterPlot, 'tabScatter', data4trajPlotNoOut, in.fname = function() return(FPDFSCATTER))
dmattek's avatar
dmattek committed
884
  
dmattek's avatar
dmattek committed
885
  ##### Hierarchical clustering
dmattek's avatar
dmattek committed
886
  callModule(clustHier, 'tabClHier', data4clust, data4trajPlotNoOut, data4stimPlot)
dmattek's avatar
dmattek committed
887
888
  
  ##### Sparse hierarchical clustering using sparcl
dmattek's avatar
dmattek committed
889
  callModule(clustHierSpar, 'tabClHierSpar', data4clust, data4trajPlotNoOut, data4stimPlot)
dmattek's avatar
dmattek committed
890

dmattek's avatar
Mod:    
dmattek committed
891
  
dmattek's avatar
dmattek committed
892
})