---
title: "meaRtools: Performing a Comprehensive Analysis of Multiple MEA Recordings"
author: "Sahar Gelfman"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{meaRtools: Performing a Comprehensive Analysis of Multiple MEA Recordings}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


meaRtools was constructed to provide core algorithms for MEA spike train analysis, feature extraction, statistical analysis and plotting of multiple MEA recordings with multiple genotypes and treatments. This vignette directs the package user on how to perform an analysis of an exemplary experiment that is made of three sequential recordings of the same plate.

## Installing the package

The meaRtools package is available for installation under the CRAN repository:

```{r install}
#install.packages( "meaRtools",repos="http://cran.us.r-project.org")
```


```{r rsetup-knitr,eval=TRUE,include=FALSE}
require(knitr)
opts_chunk$set(cache=FALSE)
```

## Loading the package

To load the package please load meaRtools and dependencies:
```{r library}
library(meaRtools)
library(plyr)
library(ggplot2)
library(reshape2)
```


## Selecting input files

Before we can perform the analysis we need to locate the recording spikeList.csv files. For this purpose we will choose three recording .csv files and a plate layout .csv file that come with the meaRtools package under /meaRtools/extdata/ :

```{r readdata}

# set path to "_spike_list.csv" files from the file path in 'filesPath'
spk_list_files<-c(system.file("extdata","exampleRecording_1012016_plate1_DIV1_spike_list.csv.gz",package = "meaRtools"),
                system.file("extdata","exampleRecording_1012016_plate1_DIV3_spike_list.csv.gz",package = "meaRtools"),
                system.file("extdata","exampleRecording_1012016_plate1_DIV4_spike_list.csv.gz",package = "meaRtools"))

# set the recording layout file "_expLog.csv"
experimental_log_file <- system.file("extdata","exampleRecording_1012016_plate1_expLog.csv.gz",package = "meaRtools")
```


## Setting output directories

Before starting the actual loading of the data, lets also set  input and output directories

```{r set-dir}
# The next command will get the directory of the csv files
data_dir<-tempdir()

# creating output under the new temp directory
print(paste0("Creating output folders under ",data_dir))

# create the output directory as /Analysis under the data_dir
output_dir<-file.path( data_dir , "Analysis" ) 
dir.create(output_dir)

# create the output directory for single recording analysis 
output_perDIV_dir<-file.path( data_dir , "Analysis/outputPerDIV" ) 
dir.create(output_perDIV_dir)

# create the output directory for R objects of analyzed recordings 
r_object_dir<-file.path( data_dir , "Analysis/R_Objects" )
dir.create(r_object_dir)

# create the output directory for log files
log.dir<-file.path( output_dir , "LogFiles" ) 
dir.create(log.dir)

# For organization sake, set a list object to hold all output directories 
analysis<-list(spikeFiles = spk_list_files, output_dir = output_dir,
           r_output_dir = r_object_dir, output_perDIV_dir = output_perDIV_dir)

```

Now let's load the recordings and create the 'spike.list' class R object

```{r spike-list}
# A loop to go over all three recording files
for (i in 1:length(spk_list_files)){
  #save title for output file name
  title<-strsplit(basename(spk_list_files[i]), ".csv")[[1]][1]
  #load plate design info for each file in the list
  plate_chem_info<-get_experimental_log_file(spk_list_files[i], experimental_log_file)
  
    # convert the spike list data to a 'spike.list' class Robject
  analysis$Robject[i]<-read_spikelist(key=title, spk_list_file=spk_list_files[i],          chem_info=plate_chem_info,r_object_dir=r_object_dir) 
}
```

## Extracting spike and burst data 
Now we have the information from each spike list file  in a new 'spike.list' R object that will be next saved in the /Analysis/Robject directory

The next step will be to construct a list of the objects for each recording and extract the features.
But first, let's load the default parameters that come with the package, each one can be set by the user. This file contains default parameters for all the functions, and we'll mention some of them here.

```{r parameters}
data("parameters")
```

You can change the timestamp of the analysis parameters so that you can track the time the analysis was done
by: parameters$timeStamp=format(Sys.time(), "%m-%d-%y_%H_%M%_%S")
For this example we are using the default which is "DATE_TIME" and will be printed in the output file names

## Extracting spike and burst features

calculate_spike_features is the first function used in the analysis pipeline, as that, this function also constructs the 'spike.list' object (called 's' here) and sets the parameters for the analysis inside the object.
We now set the defaults for an active electrode as in the parameters object, setting the minimum MFR to a lenient 1 spike in 60s and a maximum MFR of 1,000Hz. We also set the minimum of active electrodes to include a well in the analyses to 4 electrodes, which is 25% of the electrodes in a 48-well plate. 

The parameters also hold the selected algorithm for burst detection: "mi" for Maximum Interval and "si" for the Poisson Surprise algorithm in parameters$burst_type. To extract burst features we use calculate_burst_features.  For this example let's use the "ps" algorithm, so we have to set it before calling the object initializing function.

```{r spike-features}

# Construct the 'spike.list' object and calculate spike features
s<-calculate_spike_features(analysis$Robject, parameters)

# As mutual information and entropy values are considered spike features, we will calculate them here
# based on the spike data of each electrode
for (i in 1:length(s)) {
  ent_mi_data <- calculate_entropy_and_mi(s[[i]], s[[i]]$treatment, mult_factor=1.5, bin_size=0.1)
  s[[i]]$mutual_inf <- ent_mi_data[["data_dists"]][["MI"]]
  s[[i]]$entropy <- ent_mi_data [["data_dists"]][["ENT"]]
}

# Select burst algorithm
parameters$burst_type="ps"

# Detect bursts and calculate their feature statistics
s<-calculate_burst_features(s)

# Iterate through all the recordings to calculate inter-spike intervals and well level mean firing rate and add that to the 'spike.list' object

for (i in 1:length(s)) {
  s[[i]] <- calculate_isis(s[[i]])
  s[[i]]$well_stats <- compute_mean_firingrate_by_well(s[[i]])
}

```

That's it, basic spike and burst features are now stored in the 'spike.list' object (s) for each of the recordings.
We can now view them by looking into the object. For example, to view spikes for electrode B3_41 in the first recording, try the following:

```{r B3_41}
s[[1]]$spikes$B3_41
```

To view burst information for bursts calculated for electrode E7_42 of the 2nd recording, try the following:

```{r E7_42}
s[[2]]$allb$E7_42
```
beg and end stand for the sequential spike  that begins and ends the burst. IBI is the inter burst interval from the previous burst. len is the number of spikes in the burst. durn is the burst duration in seconds. mean_isis is the average inter spike intervals within this burst and SI is the surprise index, only relevant when running the poisson surprise algorithm.

#### Extracting network spike data 

To calculate network spikes we iterate over all the recordings , for each recording we  1) calculate the network spikes and 2) extract the network spike features from the spike.list object. To call network spikes, we provide the function calculate_network_spikes with several arguments: 

s[[i]] - the first is the 'spike.list' object of the recording

sur - the number of datapoints to be used in summmarizing mean network spikes (default is 100)

ns_n - the number of electrodes above which  a network spike will be called 

ns_t - the time window for calling a network spike (10ms). 

For extracting the features into the 'spike.list' object, we use IGM.summary.network.spikes and provide it with the 'spike.list' object, the calculated network spikes data, the minimum number of spikes in each electrode that we wish to consider (default is 1) and agaiun the 'sur' parameter from above.

```{r all-features}

# Iterate through all the recordings
for (i in 1:length(s)) {

  #Calculate Network Spikes
  nspikes_old <- calculate_network_spikes(s[[i]],parameters$sur, parameters$ns_n, parameters$ns_t)
  
  # Extract network spike features that will be printed later
  nspikes <- summarize_network_spikes(s[[i]],nspikes_old,ns_e = 1, parameters$sur)
  
  # Add network spike data to the 'spike.list' object
  s[[i]]$ns_all<-nspikes$ns_all
}
```

We now have all the network spike features calculated and we can look at them easily. Try running the following to see the features extracted for well B5 :

```{r en_brief}

s[[i]]$ns_all$B5$en_brief

```

#### Extracting network bursts data 

The last attribute we can extract is network bursts. To extract network bursts we do not require iterating through the 'spike.list' object, that is done automatically by calculate_network_bursts. We provide the function with several arguments alongside the 'spike.list' object: 
Sigma - the window sizes used for the analysis (10, 20 and 50ms)

min_electrodes - the minimum electrodes to call a network burst 

local_region_min_nae - to tell the algorithm if we would like to use an adaptive threshold (default is 0).

```{r nb}

   nb.list <- calculate_network_bursts(s,parameters$Sigma,
                                       parameters$min_electrodes,
                                       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(s) ){
      s[[i]]$nb_all<-nb.list$nb_all[[i]]
      s[[i]]$data.frame$nb_features<-nb.list$nb_features[[i]]
    }
    
```

## Computing spike train correlations

To compute the nature of correlated activity, we can use the spike
train tiling coefficient (STTC), proposed by Cutts and Eglen (2014, J
Neuroscience).  For example here we calculate the mean pairwise STTC,
averaged over all distinct pairs of electrodes in a well.  This is
then stored as the `mean_sttc` in the well object.

```{r sttc}
for (i in 1:length(s)) {
  s[[i]]$mean_sttc <- compute_mean_sttc_by_well(s[[i]])
}
```

For example, here are the mean correlations for all wells on the first
plate (ignoring any wells that have 0 or 1 electrode):

```{r unlist}
unlist(s[[1]]$mean_sttc)
```

For further information about the STTC measure, see the separate
vignette, `sttc.Rmd`.


## Writing and plotting single recording data  

At this point you might wonder how to produce burst feature distributions that meaRtools produces. The answer is that since the distributions derive from the burst features that were already calculated,  they are automatically constructed as part of the printing process of burst features. So, next we turn to printing the extracted features for each single recording. When printing burst features, we'll come back to producing burst feature distribution.

### printing spike data

```{r plot-summary}

# print spikes graphs (pdf format) for each recording
plot_plate_summary_for_spikes(s,analysis$output_perDIV_dir)

# write spike feature tables for each recording
suppressWarnings(write_plate_summary_for_spikes(s,analysis$output_perDIV_dir))
```

### printing burst data

As promised, the burst feature distributions are already calculated by the burst printing functions, since they are extracted from calculated burst features. The distribution features are calculated for five burst features : burst duration, IBI, nspikes (number of spikes in a burst), spikeFreq (Hz) and ISI within bursts. When running plot_plate_summary_for_bursts, the function calls calc_burst_distributions to calculate and plot those distributions for each loaded recording. The default parameters object loaded earlier, holds five objects with 7 arguments for the five distribution features: 

min_cases - the minimum number of bursts for performing the analysis

x_axis_lim - the maximum value of that feature. In this example, we perform distribution analysis for IBI. The xlimit is 20, which means that the longest IBI taken into account here would be 20s long. 

bins_in_sec - the bins in each second of IBI. Here bins_in_sec is set to 5, meaning that the IBI distribution will  be cut into 0.2s bins in a maximum of 20s. Thus, the final distribution will be made of 100 bins of 0.2s.

filter_by_min - a binary, to decide whether bursts should be filtered by a minimum value (default is 0) 

min_values - the actual minimum IBI to filter by. 

per_well - a binary argument, 1: the algorithm will group electrodes by well, and then group wells by treatment and 0 (default): electrodes will be grouped directly by treatment. 

perform - a binary argument meant to decide whether this distribution analysis should be performed at all. While   the default is to perform all five distributions, in this example we perform only the IBI distribution analysis    for all three recordings.

```{r pdfs}
# plot burst pdfs for each recording
suppressWarnings(plot_plate_summary_for_bursts(s,analysis$output_perDIV_dir,parameters))

# write burst feature tables for each recording
write_plate_summary_for_bursts(s,analysis$output_perDIV_dir)
```

### printing network spikes data

The commands below print single recording ns data. Here they are printed for the first of the three recordings

```{r ns}
i=1 
# Get plate name
basename <- strsplit(basename(s[[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,"/ns_plot.pdf"))
xyplot_network_spikes(nspikes)	
plot_active_wells_network_spikes(nspikes)
dev.off()

# write network spike data to output file
write_network_spikes_to_csv(s[[i]],nspikes,analysis$output_perDIV_dir)

# Check the graphs and csvs printed under the analysis$output_perDIV_dir path
```

## Aggregating recordings

One of the strong advantages of meaRtools is it's ability to combine the information from all the loaded recordings and use all of it when comparing between treatments. The following commands aggregate the data for spikes, bursts and network spikes. Network burst features were already aggregated automatically when we ran nb_matrix_to_feature_dfs.

```{r aggregate}
spike_features = aggregate_features(s, "spike",parameters)
ns_features = aggregate_features(s, "ns",parameters)
burst_features = aggregate_features(s, "burst",parameters)

# printing spike features nae
spike_features$nae

#Feel free to explore the spike/ns/burst and nb_features for the different features they offer
```

## Filtering inactive wells
The next step is optional, it allows the user to discard from the analysis wells that were not active in at least X% of the recordings. This percentage is also a default parameter, currently set to 50%. Thus, any well that is not considered active (at least 4 active electrodes) in at least 2 out of the 3 loaded recordings will be ignored when comparing the treatments throughout the recordings. In this example we perform the filter on spike and NB features.

```{r filtering}

# All uncalculated aEs were set previously to NA, convert all those to 0 aE before the filter
nae <- spike_features$nae
nae[is.na(nae)] <- 0

# filter spike wells
spike_features = lapply(spike_features, function(x) filter_wells( x, nae, parameters$well_min_rate, parameters$well_max_div_inactive_ratio))

# filter network burst wells
nb_features <- lapply(nb_features, function(x) filter_wells(x, nae, parameters$well_min_rate, parameters$well_max_div_inactive_ratio ))
# re-order features by well name
nb_features <- lapply(nb_features, function(x) x[order(x[,'well']),])

# printing spike features nae after filter
spike_features$nae

```

After the filter we can observe the spike features dataframe and find that well D5 was dropped from the table because it lacked activity in the first two recordings (compare to spike_features$nae before the filter was applied).

## Writing aggregated tables to files

We can now easily print all the aggregated tables of all the extracted features using one command for each attribute, These files are printed into a designated directory by the name of each activity attribute (spikes, bursts, ns, nb) under the output Analysis folder.

```{r csvs}
#write csvs 
write_features_to_files(s, spike_features, analysis$output_dir, "spikes")
write_features_to_files(s, burst_features, analysis$output_dir, "bursts")
write_features_to_files(s, ns_features, analysis$output_dir, "ns")
write_features_to_files(s, nb_features, analysis$output_dir, "nb")
```

## Testing for differences between treatments

The example recordings have three treatments (groups): treatX, treatY and untreated. to perform MW-tests, permutate the data and plot, we first need to decide which treatment we would like to compare to all the others. Here we use 'untreated' as that treatment that will be tested against treatX and treatY. However, you can also use the get_wt(s) function and it will open a tcltk window with the treatments available on the plate and will let you choose the one you're interested in, to use as argument to the testing scheme.

```{r permute}

suppressMessages(permute_features_and_plot(s, "treatX", parameters$perm_n, spike_features, "spikes", analysis$output_dir))
suppressMessages(permute_features_and_plot(s, "treatX", parameters$perm_n, burst_features, "bursts", analysis$output_dir))
suppressMessages(permute_features_and_plot(s, "treatX", parameters$perm_n, ns_features, "ns", analysis$output_dir))
suppressMessages(permute_features_and_plot(s, "treatX", parameters$perm_n, nb_features, "nb", analysis$output_dir))

```

At this point, we have extracted all features, combined the recordings, tested the differences between the treatments and printed all the results in graphs and tables. The last thing we want to perform is combining the distributions of all three recordings, testing distribution differences between treatments and printing the results. This step is easily done using one function: dist_perm. The function requires a distribution file for each burst feature, which are automatically printed by calc_burst_distributions into the /Analysis/outputPerDIV/distributionFiles folder. Aside from the distribution file, dist_perm also required the number of permutations to perform and the two treatments to be compared. dist_perm returns an objects with the following information: 
- two  distributions of the feature, one per each treatment, combined for all the recordings, the 
- two comulative distributions, as as the above
- p-values for the distribution differences (Earth Mover's Distance) and comulative distribution differences (Maximum Distance, please see meaRtools manuscript for full details)
- permutation p-values for each distribution
- original maximum distance and EMD value
The normalized distribution and it's corrsponding EMD p-value can be extracted from the dist_perm returned object and plotted as follows:
```{r dist_perm}
result <- suppressWarnings(dist_perm(paste0(analysis$output_perDIV_dir,"/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_ibi_distributions.csv"),1000,"treatX","treatY"))

plot(result$data_wt_original,col="blue",main=basename,type="l",lwd=3,xlab="IBI")
points(result$data_ko_original,col="green",type="l",lwd=3)
par(mfrow=c(1,1))  
mtext(side = 1, at = 0, line = 4,
          text = paste("P.value EMD after 1000 permutations: ",format((1-result$perm_EMD), digits = 2),sep=""),col = "black",cex= 0.9,adj=0)    
``` 

And the cumulative distribution with it's corrsponding MD p-value can be extracted as follows:
    
```{r dist_perm2}
suppressWarnings(result <- dist_perm(paste0(analysis$output_perDIV_dir,"/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_ibi_distributions.csv"),1000,"treatX","treatY"))

plot(result$data_wt,col="blue",main=basename,type="l",lwd=3,xlab="IBI")
points(result$data_ko,col="green",type="l",lwd=3)
par(mfrow=c(1,1))  
mtext(side = 1, at = 0, line = 4,
      text = paste("P.value Max distance after 1000 permutations: ",format((1-result$perm_p), digits = 3),sep=""),col = "black",cex= 0.9,adj=0)    
```    

This document provides the steps to use meaRtools's functions as a pipeline for an MEA experiment analysis lasting 3 DIVs. In the /Analysis output folder are now the full results for extracting all the features and comparing them between the treatments that we decided to test. Even with only three recordings the output files are numerous require a deep dive for better understanding the full picture of the difference between treatments. We hope you will enjoy all the capabilities that meaRtools has to offer and use it to deeply and fully understand your MEA recordings. Please read the meaRtools manuscript for detailed information about the methods and do not hesitate to contact us for any explanation that might be lacking in this document.

Goodluck!

meaRtools team

