123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325 |
- # Analysis code for MEA data
- #library(IGM.MEA)
- library(meaRtools) #attach meaRtools after IGM.MEA to get right parameter names <- fixed by getting rid of IGM.MEA, just by creating an extra function
- library(plyr)
- library(ggplot2)
- library(reshape2)
- library(gtools)
- library(pracma)
- library(svDialogs)
- CRAPels=list()
- #CODEdirectory<- folder where the code files reside
- CODEdirectory<-paste0(choose.dir(default = "", caption = "SELECT THE FOLDER WITH THE CODE"),"\\")
- #alternative, not to select this every code launch:
- ##CODEdirectory<-'Your folder with the code'. Remember / or \\ at the end
- output_folder <- choose.dir(default = "", caption = "SELECT YOUR OUTPUT FOLDER")
- options(warn=1)
- parameter_file <- paste0(CODEdirectory, 'saved.parameters')
- #Select the spike files
- my.files <- as.list(choose.files(default = "", caption = "SELECT YOUR SPIKE FILES", multi = TRUE))#, filters = Filters, index = nrow(Filters))
- my.noise <- choose.files(default = "", caption = "SELECT YOUR NOISY ELECTRODE FILE", multi = FALSE)
- my.explog <- choose.files(default = "", caption = "SELECT YOUR TREATMENT FILE (expLog)", multi = FALSE) #MULTIPLE??????
- ########## ADJUST PARAMETERS ##########
- # load parameters
- my.parameters <- readRDS(parameter_file)
- # change parameters if needed
- my.parameters$elec_min_rate <- 10/60 # minimum spikes/second for electrode to be included in analysis, e.g. 10/60
- my.parameters$elec_max_rate <- 1000 # maximum spikes/second for electrode to be included in analysis, e.g. 1000
- my.parameters$well_min_rate <- 1 # minimum number of active electrodes (here >5 spikes/minute used) for a well to be included to goodwells, e.g. 1
- my.parameters$Sigma <- c(10, 20, 50) # window sizes for network bursts analysis in ms, e.g. c(10, 20, 50)
- my.parameters$min_electrodes <- 5 # minimum number of electrodes participating in a network burst, e.g. 5
- ########## RUN ADDITIONAL FUNCTIONS ##########
- MEAlayout <- as.integer(dlgInput("Enter your MEA plate layout: 12 or 48")$res)
- if (MEAlayout == 12) {
- source(paste0(CODEdirectory, "plate_layout_Axion12_LAa.R"))
- my.array_type <- 'Axion12'
- my.pos <- paste0(CODEdirectory, 'positions_Axion_12wp_200nm_spacing.csv')
- } else if (MEAlayout == 48) {
- source(paste0(CODEdirectory, "plate_layout_Axion48_LAa.R"))
- my.array_type <- 'Axion48'
- my.pos <- paste0(CODEdirectory, 'positions_Axion_48wp_350nm_spacing.csv') #positions_Axion_48wp.csv
- } else {
- stop("Please, enter the appropriate MEA layout configuration: 12 or 48 wells per plate!")
- }
- source(paste0(CODEdirectory, "logisi_algorithm_2018_LAa.R"))
- source(paste0(CODEdirectory, "aggregate_features_algorithm_LAa.R"))
- source(paste0(CODEdirectory, "burst_summary_algorithm_LAa.R"))
- source(paste0(CODEdirectory, "github_nb_bursts_modifications_LAa.R"))
- source(paste0(CODEdirectory, "cheminfo.R"))
- ########## IMPORT ##########
- my.output <- paste0(output_folder, '/Analysis')
- dir.create(my.output)
- my.out_perDIV <- paste0(my.output, '/PerDIV')
- dir.create(my.out_perDIV)
- my.out_object <- paste0(my.output, '/R objects')
- dir.create(my.out_object)
- my.out_log <- paste0(my.output, '/LogFiles')
- dir.create(my.out_log)
- analysis <- list(spikeFiles = my.files, output_dir = my.output, r_output_dir = my.out_object, output_perDIV_dir = my.out_perDIV)
- my.exp_ID <- strsplit(basename(my.explog), split = '_')[[1]][1] #my.treatment_file
- my.div <- lapply(my.files, function(x) strsplit(basename(x), split='\\.')[[1]][1]) #my.spike_files
- my.div <- lapply(my.div, function(x) strsplit(basename(x), split = '_')[[1]][4])
- # import treatment information
- my.chem_info <- cheminfo(my.files[[1]], my.explog)
- # import noisy electrode information
- my.noise_info <- read.csv(my.noise, stringsAsFactors = FALSE)
- my.noise_div <- c()
- for (i in 1:length(my.noise_info)) {
- my.noise_day <- colnames(my.noise_info)[i]
- my.noise_div <- c(my.noise_div, my.noise_day)
- }
- # import spike files
- my.set <- list()
- for (i in 1:length(my.files)) {
- my.set[[i]] <- read_spikelist_text(my.files[[i]], my.pos, array = my.array_type, div = my.div[[i]])
- my.set[[i]]$parameters <- my.parameters
- }
- ########## ANALYSIS ##########
- # filter inactive electrodes and noisy electrodes (inactive wells not filtered)
-
- for (i in 1:length(my.set)) {
- my.wells <- my.set[[i]]$well
- my.day <- my.set[[i]]$div
- # filter inactive electrodes
- low <- which(my.set[[i]]$meanfiringrate < my.parameters$elec_min_rate)
- high <- which(my.set[[i]]$meanfiringrate > my.parameters$elec_max_rate)
- extremes <- c(low, high)
- if (!isempty(extremes)) {
- bad_ids <- names(extremes)
- CRAPels[[i]]=bad_ids
- bad_ids <- c("-", bad_ids) #MINUS IN FRONT OF THEM REMOVES NOT PRESERVES
- my.set[[i]] <- remove_spikes(my.set[[i]], bad_ids)
- }
- # filter noisy electrodes
-
- #INSERT CONDITIONAL STATEMENT HERE !! IF NO NOISY ELECTRODES----OK, PROCEED
- if (!isempty(my.set[[i]]$channels)) {
- all_ids <- names(my.set[[i]]$meanfiringrate)
- bad_ids <- c()
- if (is.element(my.day, my.noise_div)) {
- my.index <- which(my.noise_div==my.day)
- my.current <- my.noise_info[my.index]
- if (nrow(my.current)>0) {
- for (k in 1:nrow(my.noise_info)) {
- if (is.element(my.current[k,1], all_ids)) {
- bad_ids <- c(bad_ids, my.current[k,1])
- }
- }
- } else {
- message("WARNING: Seems like your noisy electrode file is empty")
- }
- }
- if (length(bad_ids) > 0) {
- bad_ids <- c("-", bad_ids)
- #my.set[[i]] <- remove_spikes(my.set[[i]], bad_ids)
- my.set[[i]] <- suppressWarnings(remove_spikes(my.set[[i]], bad_ids))
- }
- }
- my.set[[i]]$parameters <- my.parameters
- my.set[[i]]$well <- my.wells
- my.set[[i]] <- get_num_ae(my.set[[i]]) #lists for each well the # of active electrodes
- my.set[[i]]$goodwells <- names(which(my.set[[i]]$nae >= my.parameters$well_min_rate)) #writting down goodwells' labels using "names" function
- my.set[[i]]$parameters$time_stamp <- format(Sys.time(), '%m-%d-%y_%H_%M%_%S')
- my.set[[i]]$div <- my.day
- my.treatment <- my.chem_info$treatment
- my.set[[i]]$treatment <- my.treatment
- names(my.set[[i]]$treatment) <- my.chem_info$well #assignation of well labels to the tratments
- my.size <- my.chem_info$size
- my.dose <- my.chem_info$dose
- my.units <- my.chem_info$units
- my.set[[i]]$size <- as.array(my.size)
- names(my.set[[i]]$size) <- my.wells
- my.set[[i]]$dose <- as.array(my.dose)
- names(my.set[[i]]$dose) <- my.wells
- my.set[[i]]$units <- as.array(my.units)
- names(my.set[[i]]$units) <- my.wells
- }
- # remove recordings with no active electrodes
- my.inactive_rec <- c()
- my.inactive_div <- c()
- for (i in 1:length(my.set)) {
- my.index <- c()
- if (isempty(my.set[[i]]$channels)) {
- my.index <- i
- my.inactive_div <- c(my.inactive_div, my.div[[i]])
- my.inactive_rec <- c(my.inactive_rec, my.index)
- }
- }
- if (!isempty(my.inactive_rec)) {
- my.set <- my.set[-my.inactive_rec]
- my.removed_recs <- cbind((my.inactive_div), (my.inactive_rec))
- colnames(my.removed_recs) <- c('DIV','Index')
- write.csv(my.removed_recs, file = paste0(my.output,'/inactive_recordings.csv'), quote = FALSE, row.names = FALSE)
- }
- # analysis begins
- # entropy and mutual information
- for (i in 1:length(my.set)) {
- my.emi <- calculate_entropy_and_mi(my.set[[i]], my.set[[i]]$treatment, mult_factor=1.5, bin_size=0.1)
- if (isempty(my.emi[["data_dists"]][["MI"]])) {
- my.set[[i]]$mutual_inf <- NA
- } else {
- my.set[[i]]$mutual_inf <- my.emi[["data_dists"]][["MI"]]
- }
- if (isempty(my.emi[["data_dists"]][["ENT"]])) {
- my.set[[i]]$entropy <- NA
- } else {
- my.set[[i]]$entropy <- my.emi[["data_dists"]][["ENT"]]
- }
- }
- # burst features
- for (i in 1:length(my.set)) {
- current <- my.set[[i]]
- if (length(current$nspikes) > 0) {
- current$allb <-
- lapply(current$spikes, function(x) logisi.pasq.method_mod(x, cutoff = 0.1))
- current$bs <- calc_burst_summary(current, bursty_threshold = 1)
- current$bs$burst_type <- "li"
- for (j in 1:nrow(current$bs)) {
- if (current$bs[j,'nbursts'] == 0) {
- current$bs[j,'bursts_per_sec'] <- NA
- current$bs[j,'bursts_per_min'] <- NA
- current$bs[j,'mean_dur'] <- NA
- current$bs[j,'mean_spikes'] <- NA
- }
- }
- my.set[[i]] <- current
- }
- }
- # isis
- for (i in 1:length(my.set)) {
- my.set[[i]] <- calculate_isis(my.set[[i]])
- my.set[[i]]$well_stats <- compute_mean_firingrate_by_well(my.set[[i]])
- }
- # network spikes
- for (i in 1:length(my.set)) {
- nw_spikes_old <- calculate_network_spikes(my.set[[i]], my.parameters$sur, my.parameters$ns_n, my.parameters$ns_t)
- nw_spikes <- summarize_network_spikes(my.set[[i]], nw_spikes_old,ns_e = 1, my.parameters$sur)
- my.set[[i]]$ns_all <- nw_spikes$ns_all
- my.set[[i]]$nw_spikes <- nw_spikes
- }
- # networks bursts
- nb.list <- calculate_network_bursts_mod(my.set, sigmas = my.parameters$Sigma, min_electrodes = my.parameters$min_electrodes, local_region_min_nae = my.parameters$local_region_min_nAE)
- nb_features <- nb_matrix_to_feature_dfs( nb.list$nb_features_merged )
- # attach data to s object
- for (i in 1:length(my.set) ){
- my.set[[i]]$nb_all<-nb.list$nb_all[[i]]
- my.set[[i]]$data.frame$nb_features<-nb.list$nb_features[[i]]
- }
- # STTC
- for (i in 1:length(my.set)) {
- my.sttc <- c()
- my.sttc <- compute_mean_sttc_by_well(my.set[[i]])
- for (k in 1:length(my.sttc)) {
- if (is.null(my.sttc[[k]])) {
- my.sttc[[k]][1] <- NA
- }
- }
- my.set[[i]]$mean_sttc <- my.sttc
- }
- # save the data set
- saveRDS(my.set, paste0(analysis$r_output_dir, '/', my.exp_ID, format(Sys.time(), '_%m-%d-%y_%H_%M%_%S')))
- ########## EXPORT ##########
- # print spikes graphs (pdf format) for each recording
- plot_plate_summary_for_spikes(my.set,analysis$output_perDIV_dir)
- # write spike feature tables for each recording
- suppressWarnings(write_plate_summary_for_spikes(my.set,analysis$output_perDIV_dir))
- # plot burst pdfs for each recording
- suppressWarnings(plot_plate_summary_for_bursts(my.set,analysis$output_perDIV_dir,my.parameters))
- # write burst feature tables for each recording
- write_plate_summary_for_bursts_mod(my.set,analysis$output_perDIV_dir)
- # ns plot and excel commands seems weird with nspikes
- for (i in 1:length(my.set)) {
- basename <- strsplit(basename(my.set[[i]]$file), "[.]")[[1]][1]
- #Use the next commands for plotting all the ns graphs. Try opening a pdf file so that all will be printed to the same file (which is automatically done for burst features):
- pdf(file=paste0(analysis$output_perDIV_dir,"/", basename, "_ns_plot.pdf"))
- xyplot_network_spikes(my.set[[i]]$nw_spikes)
- plot_active_wells_network_spikes(my.set[[1]]$nw_spikes)
- dev.off()
- }
- # aggregate features (goodwells, emi and sttc needed)
- my.aggs <- aggregate_features(my.set, 'spike', parameters = my.parameters)
- my.aggb <- aggregate_features_mod(my.set, 'burst', parameters = my.parameters)
- my.aggns <- aggregate_features(my.set, 'ns', parameters = my.parameters)
- # export aggregated features
- write_features_to_files(my.set, my.aggs, analysis$output_dir, 'spikes')
- write_features_to_files(my.set, my.aggb, analysis$output_dir, 'bursts')
- write_features_to_files(my.set, my.aggns, analysis$output_dir, 'ns')
- write_features_to_files(my.set, nb_features, analysis$output_dir, 'nb')
- # print the following alarms if needed
- if (length(my.div) == length(my.noise_div)) {
- n <- 0
- for (i in 1:length(my.noise_div)) {
- if (is.element(my.noise_div[[i]], my.div)) {
- n <- n + 1
- }
- }
- if (n != length(my.noise_div)) {
- print(paste0((length(my.noise_div)-n), ' noise DIVs do not correspond to rec DIVs.'))
- }
- } else {
- print(paste0('Length of rec DIVs does not correspond to length of noise DIVs.'))
- }
- if (!isempty(my.inactive_rec)) {
- print(paste0('Some recordings did not include active electrodes.'))
- }
|