Comparative_MEA_dataset/Codes/meaRtools/MEA_analysis_Axion.R

# 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 .R function
library(plyr)
library(ggplot2)
library(reshape2)
library(gtools)
library(pracma)

library(svDialogs) 


#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"))

CRAPels=list()
########## 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.'))
}