server.R 37.6 KB
Newer Older
dmattek's avatar
dmattek committed
1

2

dmattek's avatar
dmattek committed
3

dmattek's avatar
dmattek committed
4
5
6
7
8
9
10
11
12
13
# This is the server logic for a Shiny web application.
# You can find out more about building applications with Shiny here:
#
# http://shiny.rstudio.com
#

library(shiny)
library(shinyjs) #http://deanattali.com/shinyjs/
library(data.table)
library(ggplot2)
dmattek's avatar
dmattek committed
14
library(gplots) # for heatmap.2
dmattek's avatar
dmattek committed
15
library(plotly)
dmattek's avatar
dmattek committed
16
17
library(d3heatmap) # for interactive heatmap
library(dendextend) # for color_branches
dmattek's avatar
Mod:    
dmattek committed
18
library(colorspace) # for palettes (ised to colour dendrogram)
dmattek's avatar
dmattek committed
19
20
21
library(RColorBrewer)
library(sparcl) # sparse hierarchical and k-means
library(scales) # for percentages on y scale
dmattek's avatar
dmattek committed
22

23
# increase file upload limit
dmattek's avatar
Added:    
dmattek committed
24
options(shiny.maxRequestSize = 80 * 1024 ^ 2)
dmattek's avatar
dmattek committed
25

dmattek's avatar
dmattek committed
26
shinyServer(function(input, output, session) {
27
  useShinyjs()
dmattek's avatar
dmattek committed
28
  
29
30
31
32
33
34
35
36
  # This is only set at session start
  # we use this as a way to determine which input was
  # clicked in the dataInBoth reactive
  counter <- reactiveValues(
    # The value of inDataGen1,2 actionButton is the number of times they were pressed
    dataGen1     = isolate(input$inDataGen1),
    dataLoadNuc  = isolate(input$inButLoadNuc)
    #dataLoadStim = isolate(input$inButLoadStim)
dmattek's avatar
dmattek committed
37
38
  )
  
dmattek's avatar
dmattek committed
39
40
41
  ####
  ## UI for side panel
  
dmattek's avatar
dmattek committed
42
  # FILE LOAD
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
  # This button will reset the inFileLoad
  observeEvent(input$inButReset, {
    reset("inFileLoadNuc")  # reset is a shinyjs function
    #reset("inButLoadStim")  # reset is a shinyjs function
  })
  
  # generate random dataset 1
  dataGen1 <- eventReactive(input$inDataGen1, {
    cat("dataGen1\n")
    
    return(userDataGen())
  })
  
  # load main data file
  dataLoadNuc <- eventReactive(input$inButLoadNuc, {
    cat("dataLoadNuc\n")
    locFilePath = input$inFileLoadNuc$datapath
    
    counter$dataLoadNuc <- input$inButLoadNuc - 1
    
    if (is.null(locFilePath) || locFilePath == '')
      return(NULL)
    else {
      return(fread(locFilePath))
    }
  })
  
dmattek's avatar
dmattek committed
70
71
72
73
74
75
76
  # This button will reset the inFileLoad
  observeEvent(input$butReset, {
    reset("inFileLoadNuc")  # reset is a shinyjs function
    #    reset("inFileStimLoad")  # reset is a shinyjs function
    
  })
  
dmattek's avatar
dmattek committed
77
78
  
  # COLUMN SELECTION
dmattek's avatar
dmattek committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
  output$varSelTrackLabel = renderUI({
    cat(file = stderr(), 'UI varSelTrackLabel\n')
    locCols = getDataNucCols()
    locColSel = locCols[locCols %like% 'rack'][1] # index 1 at the end in case more matches; select 1st
    
    selectInput(
      'inSelTrackLabel',
      'Select Track Label (e.g. objNuc_Track_ObjectsLabel):',
      locCols,
      width = '100%',
      selected = locColSel
    )
  })
  
  output$varSelTime = renderUI({
    cat(file = stderr(), 'UI varSelTime\n')
    locCols = getDataNucCols()
    locColSel = locCols[locCols %like% 'RealTime'][1] # index 1 at the end in case more matches; select 1st
    
    cat(locColSel, '\n')
    selectInput(
      'inSelTime',
      'Select time column (e.g. RealTime):',
      locCols,
      width = '100%',
      selected = locColSel
    )
  })
  
  # This is main field to select plot facet grouping
  # It's typically a column with the entire experimental description,
  # e.g. in Yannick's case it's Stim_All_Ch or Stim_All_S.
  # In Coralie's case it's a combination of 3 columns called Stimulation_...
  output$varSelGroup = renderUI({
    cat(file = stderr(), 'UI varSelGroup\n')
    locCols = getDataNucCols()
    
    if (!is.null(locCols)) {
      locColSel = locCols[locCols %like% 'ite']
      if (length(locColSel) == 0)
        locColSel = locCols[locCols %like% 'eries'][1] # index 1 at the end in case more matches; select 1st
      else if (length(locColSel) > 1) {
        locColSel = locColSel[1]
      }
      #    cat('UI varSelGroup::locColSel ', locColSel, '\n')
      selectInput(
        'inSelGroup',
        'Select one or more facet groupings (e.g. Site, Well, Channel):',
        locCols,
        width = '100%',
        selected = locColSel,
        multiple = TRUE
      )
    }
    
  })
  
  output$varSelSite = renderUI({
    cat(file = stderr(), 'UI varSelSite\n')
    
dmattek's avatar
Added:    
dmattek committed
139
140
141
142
143
144
145
146
147
148
149
150
    if (!input$chBtrackUni) {
      locCols = getDataNucCols()
      locColSel = locCols[locCols %like% 'ite'][1] # index 1 at the end in case more matches; select 1st
      
      selectInput(
        'inSelSite',
        'Select FOV (e.g. Metadata_Site or Metadata_Series):',
        locCols,
        width = '100%',
        selected = locColSel
      )
    }
dmattek's avatar
dmattek committed
151
152
153
154
155
156
157
158
159
160
  })
  
  
  
  
  output$varSelMeas1 = renderUI({
    cat(file = stderr(), 'UI varSelMeas1\n')
    locCols = getDataNucCols()
    
    if (!is.null(locCols)) {
dmattek's avatar
dmattek committed
161
162
      locColSel = locCols[locCols %like% 'objCyto_Intensity_MeanIntensity_imErkCor.*' |
                            locCols %like% 'Ratio'][1] # index 1 at the end in case more matches; select 1st
dmattek's avatar
dmattek committed
163

dmattek's avatar
dmattek committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
      selectInput(
        'inSelMeas1',
        'Select 1st measurement:',
        locCols,
        width = '100%',
        selected = locColSel
      )
    }
  })
  
  
  output$varSelMeas2 = renderUI({
    cat(file = stderr(), 'UI varSelMeas2\n')
    locCols = getDataNucCols()
    
    if (!is.null(locCols) &&
        !(input$inSelMath %in% c('', '1 / '))) {
      locColSel = locCols[locCols %like% 'objNuc_Intensity_MeanIntensity_imErkCor.*'][1] # index 1 at the end in case more matches; select 1st
dmattek's avatar
dmattek committed
182

dmattek's avatar
dmattek committed
183
184
185
186
187
188
189
190
191
192
      selectInput(
        'inSelMeas2',
        'Select 2nd measurement',
        locCols,
        width = '100%',
        selected = locColSel
      )
    }
  })
  
dmattek's avatar
dmattek committed
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
  # UI for trimming x-axis (time)
  output$uiSlTimeTrim = renderUI({
    cat(file = stderr(), 'UI uiSlTimeTrim\n')
    
    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
217
  
dmattek's avatar
dmattek committed
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
  # UI for normalization
  
  output$uiChBnorm = renderUI({
    cat(file = stderr(), 'UI uiChBnorm\n')
    
    if (input$chBnorm) {
      radioButtons(
        'rBnormMeth',
        label = 'Select method',
        choices = list('fold-change' = 'mean', 'z-score' = 'z.score')
      )
    }
  })
  
  output$uiSlNorm = renderUI({
    cat(file = stderr(), 'UI uiSlNorm\n')
    
    if (input$chBnorm) {
      locTpts  = getDataTpts()
      
      if(is.null(locTpts))
        return(NULL)
      
      locRTmin = min(locTpts)
      locRTmax = max(locTpts)
      
      sliderInput(
        'slNormRtMinMax',
        label = 'Time range for norm.',
        min = locRTmin,
        max = locRTmax,
dmattek's avatar
dmattek committed
249
250
        value = c(locRTmin, 0.1 * locRTmax), 
        step = 1
dmattek's avatar
dmattek committed
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
      )
    }
  })
  
  output$uiChBnormRobust = renderUI({
    cat(file = stderr(), 'UI uiChBnormRobust\n')
    
    if (input$chBnorm) {
      checkboxInput('chBnormRobust',
                    label = 'Robust stats',
                    FALSE)
    }
  })
  
  output$uiChBnormGroup = renderUI({
    cat(file = stderr(), 'UI uiChBnormGroup\n')
    
    if (input$chBnorm) {
      radioButtons('chBnormGroup',
dmattek's avatar
Mod:    
dmattek committed
270
271
                   label = 'Normalisation grouping',
                   choices = list('Entire dataset' = 'none', 'Per facet' = 'group', 'Per trajectory (Korean way)' = 'id'))
dmattek's avatar
dmattek committed
272
273
274
275
    }
  })
  
  
dmattek's avatar
dmattek committed
276
277
278
279
280
281
  # UI for removing outliers
  
  output$uiSlOutliers = renderUI({
    cat(file = stderr(), 'UI uiSlOutliers\n')
    
    if (input$chBoutliers) {
dmattek's avatar
Mod:    
dmattek committed
282
      
dmattek's avatar
dmattek committed
283
284
285
286
287
      sliderInput(
        'slOutliersPerc',
        label = 'Percentage of middle data',
        min = 90,
        max = 100,
dmattek's avatar
Fixed:    
dmattek committed
288
        value = 99.5, 
dmattek's avatar
dmattek committed
289
290
        step = 0.1
      )
dmattek's avatar
dmattek committed
291
      
dmattek's avatar
Mod:    
dmattek committed
292
      
dmattek's avatar
dmattek committed
293
294
295
    }
  })
  
dmattek's avatar
dmattek committed
296
297
298
299
300
301
302
303
304
  output$uiTxtOutliers = renderUI({
    if (input$chBoutliers) {
      
      p("Total tracks")
      
    }
    
  })
  
dmattek's avatar
dmattek committed
305
  
dmattek's avatar
dmattek committed
306
307
308
309
310
311
312
313
314
  ####
  ## data processing
  
  # generate random dataset 1
  dataGen1 <- eventReactive(input$inDataGen1, {
    cat("dataGen1\n")
    
    return(userDataGen())
  })
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
  
  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
    
    cat(
      "dataInBoth\ninGen1: ",
      locInGen1,
      "   prev=",
      isolate(counter$dataGen1),
      "\ninDataNuc: ",
      locInLoadNuc,
      "   prev=",
      isolate(counter$dataLoadNuc),
      # "\ninDataStim: ",
      # locInLoadStim,
      # "   prev=",
      # isolate(counter$dataLoadStim),
      "\n"
    )
    
    # isolate the checks of counter reactiveValues
    # as we set the values in this same reactive
    if (locInGen1 != isolate(counter$dataGen1)) {
      cat("dataInBoth if inDataGen1\n")
      dm = dataGen1()
      # no need to isolate updating the counter reactive values!
      counter$dataGen1 <- locInGen1
    } else if (locInLoadNuc != isolate(counter$dataLoadNuc)) {
      cat("dataInBoth if inDataLoadNuc\n")
      dm = dataLoadNuc()
      # no need to isolate updating the counter reactive values!
      counter$dataLoadNuc <- locInLoadNuc
    } else {
      cat("dataInBoth else\n")
      dm = NULL
    }
    return(dm)
  })
  
  # return column names of the main dt
dmattek's avatar
dmattek committed
364
  getDataNucCols <- reactive({
365
366
367
368
369
370
371
372
373
374
375
    cat(file = stderr(), 'getDataNucCols: in\n')
    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({
dmattek's avatar
dmattek committed
376
    cat(file = stderr(), 'dataMod\n')
377
378
    loc.dt = dataInBoth()
    
dmattek's avatar
dmattek committed
379
    if (is.null(loc.dt))
380
381
      return(NULL)
    
dmattek's avatar
Added:    
dmattek committed
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
    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'))
      {
        loc.dt[, trackObjectsLabelUni := paste(sprintf("%03d", get(input$inSelSite)),
                                               sprintf("%04d", get(input$inSelTrackLabel)),
                                               sep = "_")]
      } else if(loc.types[[input$inSelTrackLabel]] %in% c('numeric', 'integer')) {
        loc.dt[, trackObjectsLabelUni := paste(sprintf("%s", get(input$inSelSite)),
                                               sprintf("%04d", get(input$inSelTrackLabel)),
                                               sep = "_")]
      } else if(loc.types[[input$inSelSite]] %in% c('numeric', 'integer')) {
        loc.dt[, trackObjectsLabelUni := paste(sprintf("%03d", get(input$inSelSite)),
                                               sprintf("%s", get(input$inSelTrackLabel)),
                                               sep = "_")]
      } else {
        loc.dt[, trackObjectsLabelUni := paste(sprintf("%s", get(input$inSelSite)),
                                               sprintf("%s", get(input$inSelTrackLabel)),
                                               sep = "_")]
      }
dmattek's avatar
Added:    
dmattek committed
402
    } else {
dmattek's avatar
Added:    
dmattek committed
403
      loc.dt[, trackObjectsLabelUni := get(input$inSelTrackLabel)]
dmattek's avatar
Added:    
dmattek committed
404
405
    }
    
dmattek's avatar
dmattek committed
406
    
407
408
409
    return(loc.dt)
  })
  
dmattek's avatar
dmattek committed
410
411
412
413
414
  # return all unique track object labels (created in dataMod)
  # This will be used to display in UI for trajectory highlighting
  getDataTrackObjLabUni <- reactive({
    cat(file = stderr(), 'getDataTrackObjLabUni\n')
    loc.dt = dataMod()
415
    
dmattek's avatar
dmattek committed
416
417
418
419
    if (is.null(loc.dt))
      return(NULL)
    else
      return(unique(loc.dt$trackObjectsLabelUni))
420
421
  })
  
dmattek's avatar
dmattek committed
422
  # return all unique track object labels (created in dataMod)
dmattek's avatar
dmattek committed
423
  # This will be used to display in UI for trajectory highlighting
dmattek's avatar
dmattek committed
424
425
426
427
428
429
430
431
432
  getDataTrackObjLabUni_afterTrim <- reactive({
    cat(file = stderr(), 'getDataTrackObjLabUni_afterTrim\n')
    loc.dt = data4trajPlot()
    
    if (is.null(loc.dt))
      return(NULL)
    else
      return(unique(loc.dt$id))
  })
dmattek's avatar
Mod:    
dmattek committed
433
  
dmattek's avatar
dmattek committed
434
435
436
  # 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
437
438
439
  getDataTpts <- reactive({
    cat(file = stderr(), 'getDataTpts\n')
    loc.dt = dataMod()
440
    
dmattek's avatar
dmattek committed
441
442
443
444
    if (is.null(loc.dt))
      return(NULL)
    else
      return(unique(loc.dt[[input$inSelTime]]))
445
446
  })
  
dmattek's avatar
dmattek committed
447
448
449
450
451
452
453
454
455
456
457
458
459
  # return dt with cell IDs and their corresponding condition name
  # The condition is the column defined by facet groupings
  getDataCond <- reactive({
    cat(file = stderr(), 'getDataCond\n')
    loc.dt = data4trajPlot()
    
    if (is.null(loc.dt))
      return(NULL)
    else
      return(unique(loc.dt[, .(id, group)]))
    
  })
  
460
461
462
  
  # prepare data for plotting time courses
  # returns dt with these columns:
dmattek's avatar
dmattek committed
463
  #    realtime - selected from input
dmattek's avatar
dmattek committed
464
  #    y        - measurement selected from input
dmattek's avatar
dmattek committed
465
466
467
468
469
  #               (can be a single column or result of an operation on two cols)
  #    id       - trackObjectsLabelUni (created in dataMod)
  #    group    - grouping variable for facetting from input
  #    mid.in   - column with trajectory selection status from the input file or
  #               highlight status from UI
470
  data4trajPlot <- reactive({
dmattek's avatar
dmattek committed
471
    cat(file = stderr(), 'data4trajPlot\n')
472
473
    
    loc.dt = dataMod()
dmattek's avatar
dmattek committed
474
    if (is.null(loc.dt))
475
476
477
      return(NULL)
    
    
dmattek's avatar
dmattek committed
478
    if (input$inSelMath == '')
479
480
481
482
483
484
485
486
487
      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)
    
    # create expression for parsing
    # creates a merged column based on other columns from input
    # used for grouping of plot facets
dmattek's avatar
dmattek committed
488
489
490
491
    if(length(input$inSelGroup) == 0)
      return(NULL)
    loc.s.gr = sprintf("paste(%s, sep=';')",
                       paste(input$inSelGroup, sep = '', collapse = ','))
492
493
494
    
    loc.s.rt = input$inSelTime
    
dmattek's avatar
dmattek committed
495
496
497
498
    # Assign tracks selected for highlighting in UI
    loc.tracks.highlight = input$inSelHighlight
    locBut = input$chBhighlightTraj
    
dmattek's avatar
Added:    
dmattek committed
499
500
    
    # Find column names with position
dmattek's avatar
Mod:    
dmattek committed
501
502
    loc.s.pos.x = names(loc.dt)[names(loc.dt) %like% c('.*ocation.*X') | names(loc.dt) %like% c('.*os.x')]
    loc.s.pos.y = names(loc.dt)[names(loc.dt) %like% c('.*ocation.*Y') | names(loc.dt) %like% c('.*os.y')]
dmattek's avatar
Added:    
dmattek committed
503
504
505
506
507
508
    
    if (length(loc.s.pos.x) == 1 & length(loc.s.pos.y) == 1)
      locPos = TRUE
    else
      locPos = FALSE
    
509
510
511
    # if dataset contains column mid.in with trajectory filtering status,
    # then, include it in plotting
    if (sum(names(loc.dt) %in% 'mid.in') > 0) {
dmattek's avatar
Added:    
dmattek committed
512
513
514
515
516
517
      if (locPos) # position columns present
        loc.out = loc.dt[, .(
          y = eval(parse(text = loc.s.y)),
          id = trackObjectsLabelUni,
          group = eval(parse(text = loc.s.gr)),
          realtime = eval(parse(text = loc.s.rt)),
dmattek's avatar
Mod:    
dmattek committed
518
519
          pos.x = get(loc.s.pos.x),
          pos.y = get(loc.s.pos.y),
dmattek's avatar
Added:    
dmattek committed
520
          mid.in = mid.in
dmattek's avatar
Mod:    
dmattek committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
        )] else
          loc.out = loc.dt[, .(
            y = eval(parse(text = loc.s.y)),
            id = trackObjectsLabelUni,
            group = eval(parse(text = loc.s.gr)),
            realtime = eval(parse(text = loc.s.rt)),
            mid.in = mid.in
          )]
        
        
        
        
        # add 3rd level with status of track selection
        # to a column with trajectory filtering status
        if (locBut) {
          loc.out[, mid.in := ifelse(id %in% loc.tracks.highlight, 'SELECTED', mid.in)]
        }
        
539
    } else {
dmattek's avatar
Added:    
dmattek committed
540
541
542
543
544
      if (locPos) # position columns present
        loc.out = loc.dt[, .(
          y = eval(parse(text = loc.s.y)),
          id = trackObjectsLabelUni,
          group = eval(parse(text = loc.s.gr)),
dmattek's avatar
Mod:    
dmattek committed
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
          realtime = eval(parse(text = loc.s.rt)),
          pos.x = get(loc.s.pos.x),
          pos.y = get(loc.s.pos.y)
        )] else
          loc.out = loc.dt[, .(
            y = eval(parse(text = loc.s.y)),
            id = trackObjectsLabelUni,
            group = eval(parse(text = loc.s.gr)),
            realtime = eval(parse(text = loc.s.rt))
          )]
        
        
        # add a column with status of track selection
        if (locBut) {
          loc.out[, mid.in := ifelse(id %in% loc.tracks.highlight, 'SELECTED', 'NOT SEL')]
        }
561
    }
562
    
dmattek's avatar
Added:    
dmattek committed
563
564
    # add XY location if present in the dataset
    
dmattek's avatar
dmattek committed
565
566
    # remove NAs
    loc.out = loc.out[complete.cases(loc.out)]
dmattek's avatar
Mod:    
dmattek committed
567
    
dmattek's avatar
dmattek committed
568
569
570
571
    # Trim x-axis (time)
    if(input$chBtimeTrim) {
      loc.out = loc.out[realtime >= input$slTimeTrim[[1]] & realtime <= input$slTimeTrim[[2]] ]
    }
dmattek's avatar
dmattek committed
572
573
    
    # Normalization
dmattek's avatar
dmattek committed
574
    # F-n myNorm adds additional column with .norm suffix
dmattek's avatar
dmattek committed
575
576
577
578
579
580
581
582
583
584
585
586
    if (input$chBnorm) {
      loc.out = myNorm(
        in.dt = loc.out,
        in.meas.col = 'y',
        in.rt.col = 'realtime',
        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
587
588
      # Column with normalized data is renamed to the original name
      # Further code assumes column name y produced by data4trajPlot
dmattek's avatar
dmattek committed
589
590
591
592
      loc.out[, y := NULL]
      setnames(loc.out, 'y.norm', 'y')
    }
    
dmattek's avatar
dmattek committed
593
594
595
596
597
598
    ##### MOD HERE
    ## display number of filtered tracks in textUI: uiTxtOutliers
    ## How? 
    ## 1. through reactive values?
    ## 2. through additional comumn to tag outliers?
    
dmattek's avatar
dmattek committed
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
    # Remove outliers
    # 1. Scale all points (independently per track)
    # 2. Pick time points that exceed the bounds
    # 3. Identify IDs of outliers
    # 4. Select cells that don't have these IDs
    
    cat('Ncells orig = ', length(unique(loc.out$id)), '\n')
    
    if (input$chBoutliers) {
      loc.out[, y.sc := scale(y)]  
      loc.tmp = loc.out[ y.sc < quantile(y.sc, (1 - input$slOutliersPerc * 0.01)*0.5) | 
                           y.sc > quantile(y.sc, 1 - (1 - input$slOutliersPerc * 0.01)*0.5)]
      loc.out = loc.out[!(id %in% unique(loc.tmp$id))]
      loc.out[, y.sc := NULL]
    }
    
    cat('Ncells trim = ', length(unique(loc.out$id)), '\n')
dmattek's avatar
Mod:    
dmattek committed
616
    
dmattek's avatar
dmattek committed
617
    return(loc.out)
dmattek's avatar
dmattek committed
618
619
  })
  
dmattek's avatar
dmattek committed
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
  
  
  # prepare data for clustering
  # return a matrix with:
  # cells as columns
  # time points as rows
  data4clust <- reactive({
    cat(file = stderr(), 'data4clust\n')
    
    loc.dt = data4trajPlot()
    if (is.null(loc.dt))
      return(NULL)
    
    loc.out = dcast(loc.dt, id ~ realtime, value.var = 'y')
    loc.rownames = loc.out$id
    
dmattek's avatar
Mod:    
dmattek committed
636
    
dmattek's avatar
dmattek committed
637
638
639
    loc.out = as.matrix(loc.out[, -1])
    rownames(loc.out) = loc.rownames
    return(loc.out)
dmattek's avatar
Mod:    
dmattek committed
640
  }) 
dmattek's avatar
dmattek committed
641
  
dmattek's avatar
dmattek committed
642
643
  
  # get cell IDs with cluster assignments depending on dendrogram cut
dmattek's avatar
dmattek committed
644
645
  getDataCl = function(in.dend, in.k, in.ids) {
    cat(file = stderr(), 'getDataCl \n')
dmattek's avatar
Mod:    
dmattek committed
646
    
dmattek's avatar
dmattek committed
647
    loc.dt.cl = data.table(id = in.ids,
dmattek's avatar
dmattek committed
648
649
650
                           cl = cutree(as.dendrogram(in.dend), k = in.k))
  }
  
dmattek's avatar
dmattek committed
651
652
  ####
  ## UI for trajectory plot
dmattek's avatar
dmattek committed
653
  
dmattek's avatar
dmattek committed
654
655
  output$varSelHighlight = renderUI({
    cat(file = stderr(), 'UI varSelHighlight\n')
dmattek's avatar
dmattek committed
656
    
dmattek's avatar
dmattek committed
657
658
659
    locBut = input$chBhighlightTraj
    if (!locBut)
      return(NULL)
dmattek's avatar
dmattek committed
660
    
dmattek's avatar
dmattek committed
661
    loc.v = getDataTrackObjLabUni()
dmattek's avatar
dmattek committed
662
    if (!is.null(loc.v)) {
663
      selectInput(
dmattek's avatar
dmattek committed
664
665
666
        'inSelHighlight',
        'Select one or more rajectories:',
        loc.v,
667
        width = '100%',
dmattek's avatar
dmattek committed
668
        multiple = TRUE
669
      )
dmattek's avatar
dmattek committed
670
671
672
    }
  })
  
dmattek's avatar
Mod:    
dmattek committed
673
  callModule(modTrajPlot, 'modTrajPlot', data4trajPlot)
dmattek's avatar
dmattek committed
674
  
dmattek's avatar
Added:    
dmattek committed
675
676
  ###### Box-plot
  callModule(tabBoxPlot, 'tabBoxPlot', data4trajPlot)
dmattek's avatar
dmattek committed
677
  
dmattek's avatar
dmattek committed
678
679
  
  
dmattek's avatar
dmattek committed
680
681
682
  ###### Scatter plot
  callModule(tabScatterPlot, 'tabScatter', data4trajPlot)
  
dmattek's avatar
dmattek committed
683
684
685
686
687
688
689
690
691
692
693
  ##### Hierarchical clustering
  
  output$uiPlotHierClSel = renderUI({
    if(input$chBPlotHierClSel) {
      selectInput('inPlotHierClSel', 'Select clusters to display', 
                  choices = seq(1, input$inPlotHierNclust, 1),
                  multiple = TRUE, 
                  selected = 1)
    }
  })
  
dmattek's avatar
Mod:    
dmattek committed
694
695
  # perform hierarchical clustering and return dendrogram coloured according to cutree
  # branch coloring performed using dendextend package
dmattek's avatar
dmattek committed
696
  userFitDendHier <- reactive({
dmattek's avatar
Mod:    
dmattek committed
697
698
    cat(file = stderr(), 'userFitDendHier \n')

dmattek's avatar
dmattek committed
699
700
    dm.t = data4clust()
    if (is.null(dm.t)) {
dmattek's avatar
Mod:    
dmattek committed
701
      return(NULL)
dmattek's avatar
dmattek committed
702
703
704
705
706
707
708
    }
    
    cl.dist = dist(dm.t, method = s.cl.diss[as.numeric(input$selectPlotHierDiss)])
    cl.hc = hclust(cl.dist, method = s.cl.linkage[as.numeric(input$selectPlotHierLinkage)])
    cl.lev = rev(row.names(dm.t))
    
    dend <- as.dendrogram(cl.hc)
dmattek's avatar
Mod:    
dmattek committed
709
710
711
    dend <- color_branches(dend, 
                           col = rainbow_hcl, # make sure that n here equals max in the input$inPlotHierNclust slider
                           k = input$inPlotHierNclust)
dmattek's avatar
dmattek committed
712
713
714
715
    
    return(dend)
  })
  
dmattek's avatar
Mod:    
dmattek committed
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
  # prepares a table with cluster numbers in 1st column and colour assignments in 2nd column
  # the number of rows is determined by dendrogram cut
  getClCol <- function(in.dend, in.k) {
    
    loc.col_labels <- get_leaves_branches_col(in.dend)
    loc.col_labels <- loc.col_labels[order(order.dendrogram(in.dend))]
    
    return(unique(
      data.table(cl.no = dendextend::cutree(in.dend, k = in.k, order_clusters_as_data = TRUE),
                 cl.col = loc.col_labels)))
  }
  
  # returns table prepared with f-n getClCol
  # for hierarchical clustering
  getClColHier <- reactive({
    cat(file = stderr(), 'getClColHier \n')
    
    loc.dend = userFitDendHier()
    if (is.null(loc.dend))
      return(NULL)
    
    return(getClCol(loc.dend, input$inPlotHierNclust))
  })
  

  
    # Function instead of reactive as per:
dmattek's avatar
dmattek committed
743
744
745
746
747
748
749
750
751
752
753
754
  # http://stackoverflow.com/questions/26764481/downloading-png-from-shiny-r
  # This function is used to plot and to downoad a pdf
  plotHier <- function() {
    
    loc.dm = data4clust()
    if (is.null(loc.dm))
      return(NULL)
    
    loc.dend <- userFitDendHier()
    if (is.null(loc.dend))
      return(NULL)
    
dmattek's avatar
Mod:    
dmattek committed
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
    loc.p = myPlotHeatmap(loc.dm,
                  loc.dend, 
                  palette.arg = input$selectPlotHierPalette, 
                  palette.rev.arg = input$inPlotHierRevPalette, 
                  dend.show.arg = input$selectPlotHierDend, 
                  key.show.arg = input$selectPlotHierKey, 
                  margin.x.arg = input$inPlotHierMarginX, 
                  margin.y.arg = input$inPlotHierMarginY, 
                  nacol.arg = input$inPlotHierNAcolor, 
                  font.row.arg = input$inPlotHierFontX, 
                  font.col.arg = input$inPlotHierFontY, 
                  title.arg = paste(
                    "Distance measure: ",
                    s.cl.diss[as.numeric(input$selectPlotHierDiss)],
                    "\nLinkage method: ",
                    s.cl.linkage[as.numeric(input$selectPlotHierLinkage)]
                  ))
dmattek's avatar
dmattek committed
772
773
774
775
776
    
    return(loc.p)
  }
  
  
dmattek's avatar
Mod:    
dmattek committed
777
778
779
780
781
  # prepare data for plotting trajectories per cluster
  # outputs dt as data4trajPlot but with an additional column 'cl' that holds cluster numbers
  # additionally some clusters are omitted according to manual selection
  data4trajPlotCl <- reactive({
    cat(file = stderr(), 'data4trajPlotCl: in\n')
dmattek's avatar
dmattek committed
782
    
dmattek's avatar
Mod:    
dmattek committed
783
    loc.dt = data4trajPlot()
dmattek's avatar
dmattek committed
784
785
    
    if (is.null(loc.dt)) {
dmattek's avatar
Mod:    
dmattek committed
786
      cat(file = stderr(), 'data4trajPlotCl: dt is NULL\n')
dmattek's avatar
dmattek committed
787
788
789
      return(NULL)
    }
    
dmattek's avatar
Mod:    
dmattek committed
790
    cat(file = stderr(), 'data4trajPlotCl: dt not NULL\n')
dmattek's avatar
dmattek committed
791
792
    
    # get cellIDs with cluster assignments based on dendrogram cut
dmattek's avatar
dmattek committed
793
    loc.dt.cl = getDataCl(userFitDendHier(), input$inPlotHierNclust, getDataTrackObjLabUni_afterTrim())
dmattek's avatar
dmattek committed
794
795
796
    loc.dt = merge(loc.dt, loc.dt.cl, by = 'id')
    
    # display only selected clusters
dmattek's avatar
Mod:    
dmattek committed
797
798
    if(input$chBPlotHierClSel)
      loc.dt = loc.dt[cl %in% input$inPlotHierClSel]
dmattek's avatar
dmattek committed
799
    
dmattek's avatar
Mod:    
dmattek committed
800
801
    return(loc.dt)    
  })
dmattek's avatar
dmattek committed
802
  
dmattek's avatar
Mod:    
dmattek committed
803
804
805
806
807
  callModule(modTrajPlot, 'modPlotHierTraj', 
             in.data = data4trajPlotCl, 
             in.facet = 'cl',  
             in.facet.color = getClColHier,
             in.fname = paste0('clust_hierch_tCourses_',
dmattek's avatar
Mod:    
dmattek committed
808
809
810
                                                                            s.cl.diss[as.numeric(input$selectPlotHierDiss)],
                                                                            '_',
                                                                            s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.pdf'))
dmattek's avatar
dmattek committed
811
  
dmattek's avatar
Added:    
dmattek committed
812
  # download a list of cellIDs with cluster assignments
dmattek's avatar
dmattek committed
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
  output$downCellCl <- downloadHandler(
    filename = function() {
      paste0('clust_hierch_data_',
             s.cl.diss[as.numeric(input$selectPlotHierDiss)],
             '_',
             s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.csv')
    },
    
    content = function(file) {
      write.csv(x = getDataCl(userFitDendHier(), input$inPlotHierNclust, getDataTrackObjLabUni_afterTrim()), file = file, row.names = FALSE)
    }
  )
  
  output$downCellClSpar <- downloadHandler(
    filename = function() {
      paste0('clust_hierchSpar_data_',
             s.cl.spar.diss[as.numeric(input$selectPlotHierSparDiss)],
             '_',
             s.cl.spar.linkage[as.numeric(input$selectPlotHierSparLinkage)], '.csv')
    },
    
    content = function(file) {
      write.csv(x = getDataCl(userFitDendHierSpar(), input$inPlotHierSparNclust, getDataTrackObjLabUni_afterTrim()), file = file, row.names = FALSE)
    }
  )
  
  
dmattek's avatar
Mod:    
dmattek committed
840
841
842
843
844
845
  # callModule(downCellCl, 'downDataHier', paste0('clust_hierch_data_',
  #                                               s.cl.diss[as.numeric(input$selectPlotHierDiss)],
  #                                               '_',
  #                                               s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.csv'),
  #            getDataCl(userFitDendHier, input$inPlotHierNclust, getDataTrackObjLabUni_afterTrim))
  # 
dmattek's avatar
Added:    
dmattek committed
846
  
dmattek's avatar
Mod:    
dmattek committed
847
848
849
850
851
852
853
854
855
856
857
  output$downloadDataClean <- downloadHandler(
    filename = 'tCoursesSelected_clean.csv',
    content = function(file) {
      write.csv(data4trajPlot(), file, row.names = FALSE)
    }
  )
  
  
  # prepare data for barplot with distribution of items per condition  
  data4clDistPlot <- reactive({
    cat(file = stderr(), 'data4clDistPlot: in\n')
dmattek's avatar
dmattek committed
858
859
    
    # get cell IDs with cluster assignments depending on dendrogram cut
dmattek's avatar
Mod:    
dmattek committed
860
    loc.dend <- userFitDendHier()
dmattek's avatar
dmattek committed
861
862
863
864
865
866
867
868
869
    if (is.null(loc.dend)) {
      cat(file = stderr(), 'plotClDist: loc.dend is NULL\n')
      return(NULL)
    }
    
    loc.dt.cl = data.table(id = getDataTrackObjLabUni_afterTrim(),
                           cl = cutree(as.dendrogram(loc.dend), k = input$inPlotHierNclust))
    
    
dmattek's avatar
dmattek committed
870
    # get cellIDs with condition name
dmattek's avatar
Mod:    
dmattek committed
871
    loc.dt.gr = getDataCond()
dmattek's avatar
dmattek committed
872
873
874
875
876
877
878
879
    if (is.null(loc.dt.gr)) {
      cat(file = stderr(), 'plotClDist: loc.dt.gr is NULL\n')
      return(NULL)
    }
    
    loc.dt = merge(loc.dt.cl, loc.dt.gr, by = 'id')
    
    # display only selected clusters
dmattek's avatar
Mod:    
dmattek committed
880
881
    if(input$chBPlotHierClSel)
      loc.dt = loc.dt[cl %in% input$inPlotHierClSel]
dmattek's avatar
dmattek committed
882
883
884
    
    loc.dt.aggr = loc.dt[, .(nCells = .N), by = .(group, cl)]
    
dmattek's avatar
Mod:    
dmattek committed
885
    return(loc.dt.aggr)
dmattek's avatar
dmattek committed
886
    
dmattek's avatar
Mod:    
dmattek committed
887
  })
dmattek's avatar
dmattek committed
888
  
dmattek's avatar
Mod:    
dmattek committed
889

dmattek's avatar
dmattek committed
890
891
892
893
894
895
896
897
898
899
900
901
902
  #  Hierarchical - display heatmap
  getPlotHierHeatMapHeight <- function() {
    return (input$inPlotHierHeatMapHeight)
  }
  
  output$outPlotHier <- renderPlot({
    locBut = input$butPlotHierHeatMap
    
    if (locBut == 0) {
      cat(file = stderr(), 'outPlotHier: Go button not pressed\n')
      
      return(NULL)
    }
dmattek's avatar
Mod:    
dmattek committed
903

dmattek's avatar
dmattek committed
904
905
906
907
    plotHier()
  }, height = getPlotHierHeatMapHeight)
  
  
dmattek's avatar
dmattek committed
908
909
910
911
  #  Hierarchical - Heat Map - download pdf
  callModule(downPlot, "downPlotHier",       paste0('clust_hierch_heatMap_',
                                                    s.cl.diss[as.numeric(input$selectPlotHierDiss)],
                                                    '_',
dmattek's avatar
Mod:    
dmattek committed
912
                                                    s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.png'), plotHier)
dmattek's avatar
dmattek committed
913

dmattek's avatar
Mod:    
dmattek committed
914
915
916
917
  callModule(modClDistPlot, 'hierClDistPlot', 
             in.data = data4clDistPlot,
             in.cols = getClColHier,
             in.fname = paste0('clust_hierch_clDist_',
dmattek's avatar
Mod:    
dmattek committed
918
919
920
921
                    s.cl.diss[as.numeric(input$selectPlotHierDiss)],
                    '_',
                    s.cl.linkage[as.numeric(input$selectPlotHierLinkage)], '.pdf'))
  
dmattek's avatar
dmattek committed
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
  
  ##### Sparse hierarchical clustering using sparcl
  
  # UI for advanced options
  output$uiPlotHierSparNperms = renderUI({
    if (input$inHierSparAdv)
      sliderInput(
        'inPlotHierSparNperms',
        'Number of permutations',
        min = 1,
        max = 20,
        value = 1,
        step = 1,
        ticks = TRUE
      )
  })
  
  # UI for advanced options
  output$uiPlotHierSparNiter = renderUI({
    if (input$inHierSparAdv)
      sliderInput(
        'inPlotHierSparNiter',
        'Number of iterations',
        min = 1,
        max = 50,
        value = 1,
        step = 1,
        ticks = TRUE
      )
  })
  
  output$uiPlotHierSparClSel = renderUI({
    if(input$chBPlotHierSparClSel) {
      selectInput('inPlotHierSparClSel', 'Select clusters to display', 
                  choices = seq(1, input$inPlotHierSparNclust, 1),
                  multiple = TRUE, 
                  selected = 1)
    }
  })
  
dmattek's avatar
Mod:    
dmattek committed
962
  
dmattek's avatar
dmattek committed
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
  getPlotHierSparHeatMapHeight <- function() {
    return (input$inPlotHierSparHeatMapHeight)
  }
  
  userFitHierSpar <- reactive({
    dm.t = data4clust()
    if (is.null(dm.t)) {
      return()
    }
    
    #cat('rownames: ', rownames(dm.t), '\n')
    
    perm.out <- HierarchicalSparseCluster.permute(
      dm.t,
      wbounds = NULL,
      nperms = ifelse(input$inHierSparAdv, input$inPlotHierSparNperms, 1),
      dissimilarity = s.cl.spar.diss[as.numeric(input$selectPlotHierSparDiss)]
    )
    
    sparsehc <- HierarchicalSparseCluster(
      dists = perm.out$dists,
      wbound = perm.out$bestw,
      niter = ifelse(input$inHierSparAdv, input$inPlotHierSparNiter, 1),
      method = s.cl.spar.linkage[as.numeric(input$selectPlotHierSparLinkage)],
      dissimilarity = s.cl.spar.diss[as.numeric(input$selectPlotHierSparDiss)]
    )
    return(sparsehc)
  })
  
  
  userFitDendHierSpar <- reactive({
    sparsehc = userFitHierSpar()
    if (is.null(sparsehc)) {
      return()
    }
    
    dend <- as.dendrogram(sparsehc$hc)
dmattek's avatar
Mod:    
dmattek committed
1000
    dend <- color_branches(dend,