# Laura Aarnos 190515
# Plotting burst detection results

# attach packages from library

#library(IGM.MEA)
library(Hmisc)
library(meaRtools) #attach meaRtools after IGM.MEA to get right parameter names

# import data

# location of RDS file (e.g. 'Z:\Public user documents\User\Folders between\Analysis\R objects\Filename')
my.input <- 'Z:/Public user documents/paavitan/Varaosat2/Lundbeck erilaistus/MEA Analysis/Analysed files/meaRtools/Lundbeck_meaRtools/LundbeckE6/N2579_final/Analysis/R objects/N2579_02-26-19_14_17_31'
# import function
burst.set <- readRDS(my.input)

# instructions to derive data of a single electrode:
# with a click, open burst.set located in the environment
# see which [[list number]] corresponds to desired DIV
# set the information below

my.e <- burst.set[[20]]$allb$C3_14 # burst.set[[i]]$allb$e, where i is the list number, e is electrode name
my.bs <- burst.set[[20]]$bs['C3_14',] # burst.set[[i]]$bs['e',]
my.spikes <- burst.set[[20]]$spikes$C3_14 # burst.set[[i]]$spikes$e
name <- 'DIV 120, C3_14' # 'name', name is used in plots
burst.set[[20]]$rec_time # bursts.set[[i]]$rec_time, printing recording time start and end in seconds (rec time is based on first and last spike detected)

# printing mean and median burst durations

mean(my.e[,'durn']) # printing mean burst duration
median(my.e[,'durn']) # printing median burst duration

# plot boxplots of IBI, number of spikes in a burst, burst duration and mean ISIs in burst

par(mfrow = c(2,2)) # setting number of plots in display
boxplot(my.e[,'ibi'],
        names = c(' '), main = paste0(name, ', IBI'), ylab = 'IBI (s)',
        col = c('lightseagreen'))
boxplot(my.e[,'len'],
        names = c(' '), main = paste0(name, ', spikes in burst'), ylab = 'Number of spikes in burst',
        col = c('lightseagreen'))
boxplot(my.e[,'durn'],
        names = c(' '), main = paste0(name, ', burst duration'), ylab = 'Burst duration (s)',
        col = c('lightseagreen'))
boxplot(my.e[,'mean_isis'],
        names = c(' '), main = paste0(name, ', mean ISIs'), ylab = 'Mean ISIs (s)',
        col = c('lightseagreen'))

# plot histogram: % of in-burst spikes of all spikes, number of bursts

par(mfrow = c(1,2)) # setting number of plots in display
my.nbursts <- my.bs[,'nbursts']
barplot(my.nbursts,
        main = name, xlab = 'Number of bursts', ylim = c(0,max(my.nbursts)+10),
        col = c('lightseagreen'), beside = TRUE)
my.per_spikes <- my.bs[,'per_spikes_in_burst']
barplot(my.per_spikes,
        main = name, xlab = '% of spikes within bursts', ylim = c(0,100),
        col = c('lightseagreen'), beside = TRUE)

# plotting burst detection - set start and end time below

my.y <- c(rep(1, length(my.spikes))) # creating ready-to-plot spikes from time stamps
my.logISI <- cbind(my.e[,'beg'], my.e[,'end']) # retrieving burst timing info

par(mfrow = c(1,1)) # setting number of plots in display

my.start <- 0 # SET START TIME FOR PLOT
my.end <- 200 # SET END TIME FOR PLOT
my.lwd <- 3 # line width for burst line
my.nx = 5 # number of minor tick to draw between major ticks, good number depends on your plotted time interval

# plot spikes of the specified time interval
plot(my.spikes, my.y, xlim = c(my.start, my.end), ylim = c(0,2.5), type = 'h', xlab = 'Spike times (s)', ylab = '', main = name, yaxt = 'n')
minor.tick(nx = my.nx, ny = 0, tick.ratio = 0.7)

#plot bursts of the specified time interval
for (i in 1:nrow(my.logISI)) {
  lines(my.spikes[my.logISI[i,]],c(1.1,1.1), col = 'lightseagreen', lwd = my.lwd)
}

# if you wish to save your plot, click export on top of the plot