
##Install packages
install.packages("scbursts")
install.packages("MASS")
install.packages("extremevalues")

##Set working domain and load packages
setwd("/Users/christiantessier/Desktop/WD/Figure4")
library(scbursts)
packageVersion("scbursts")
library(MASS)
library(extremevalues)

##Import file, convert to dwells and define bursts
infile <- "Figure4_QX10uM_fileC.evt"
transitions <- evt.read(infile)
dwells <- evt.to_dwells(transitions)
dwells_c <- risetime.correct_gaussian(Tr=35.0052278, dwells, unit="us")
bursts <- bursts.defined_by_tcrit(dwells_c, 0.004462, units="s")

##Remove bursts with less than 3 events
threeplus<-function(seg){return(segment.count_open(seg)>=3)}
bursts2<-bursts.select(bursts,threeplus)
bursts2combined <- bursts.space_out(bursts2, sep_factor = 1000)
bursts2t <- bursts.defined_by_tcrit(bursts2combined, 0.004462, units="s")
bursts2s <- bursts.select(bursts2t, segment.verify)
burst2popen<-bursts.popens(bursts2s)

##Remove outlier popens with extremevalues
val<-getOutliers(burst2popen,method="I")
outliers <-function(seg){segment.popen(seg)<val$limit[2]&&segment.popen(seg)>val$limit[1]}
bursts3 <-bursts.select(bursts2s,outliers)
bursts3combined <- bursts.space_out(bursts3, sep_factor = 1000)
bursts3t <- bursts.defined_by_tcrit(bursts3combined, 0.004462, units="s")
bursts3s <- bursts.select(bursts3t, segment.verify)
burst3popen <- bursts.popens(bursts3s)

##Remove data beyond 2SD of the mean
fit<-fitdistr(burst3popen, "normal")
para<-fit$estimate
upper_bound<-para[1] + (2*para[2])
lower_bound<-para[1] - (2*para[2])
removebeyond2sd<-function(seg){segment.popen(seg)<upper_bound&&segment.popen(seg)>lower_bound}
bursts4 <-bursts.select(bursts3s,removebeyond2sd)
bursts4combined <- bursts.space_out(bursts4, sep_factor = 1000)
bursts4t <- bursts.defined_by_tcrit(bursts4combined, 0.004462, units="s")
bursts4s <- bursts.select(bursts4t, segment.verify)
burst10popen <- bursts.popens(bursts4)

##Write Bursts to DWT files
dwt.write(bursts4s, file="Figure4_QX10uM_fileC_Analyzed.dwt")
