IGM.MEA: Performing a Comprehensive Analysis of Multiple MEA Recordings

IGM.MEA was constructed to perform 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 IGM.MEA package is available for installation under the CRAN repository:

#install.packages( "IGM.MEA",repos="http://cran.us.r-project.org")

Loading the package

To load the package please load IGM.MEA and dependencies:


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 IGM.MEA package under /IGM.MEA/extdata/ :

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

# set the recording layout file "_expLog.csv"
ExperimentalLogFile <- system.file("extdata","exampleRecording_1012016_plate1_expLog.csv.gz",package = "IGM.MEA")

Setting output directories

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

# The next command will get the directory of the csv files

# create the output directory as /Analysis under the data.dir
output.dir<-paste0( data.dir , "/Analysis" ) 
suppressWarnings(  dir.create(output.dir) )

# create the output directory for single recording analysis 
output.perDIV.dir<-paste0( data.dir , "/Analysis/outputPerDIV" ) 
suppressWarnings(  dir.create(output.perDIV.dir) )

# create the output directory for R objects of analyzed recordings 
Robject.dir<-paste0( data.dir , "/Analysis/R_Objects" )
suppressWarnings(  dir.create(Robject.dir) )

# create the output directory for log files
log.dir<-paste0( output.dir , "/LogFiles" ) 
suppressWarnings(  dir.create(log.dir) )

# For organization sake, set a list object to hold all output directories 
analysis<-list(spikeFiles = spkListFiles, output.dir = output.dir, Routput.dir = Robject.dir, output.perDIV.dir = output.perDIV.dir)

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

# A loop to go over all three recording files
for (i in 1:length(spkListFiles)){
  #save title for output file name
  title<-strsplit(basename(spkListFiles[i]), ".csv")[[1]][1]
  #load plate design info for each file in the list
  plate.chem.info<-get.experimental.log.file(spkListFiles[i], ExperimentalLogFile)

    # convert the spike list data to a 'spike.list' class Robject
  analysis$Robject[i]<-read.spikelist(key=title, spkListFile=spkListFiles[i],    chem.info=plate.chem.info,Robject.dir=Robject.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.


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.

# Select burst algorithm

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

# Detect bursts and calculate their feature statistics

# 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 <- IGM.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:


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


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.

# 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

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 :


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).

   nb.list <- calculate.network.bursts(s,parameters$Sigma,

    nb.features <- NB.matrix.to.feature.dfs( nb.list$nb.features.merged )

    # attach data to s object
    for (i in 1:length(s) ){

Writing and plotting single recording data

At this point you might wonder how to produce burst feature distributions that IGM.MEA 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

# print spikes graphs (pdf format) for each recording

# write spike feature tables for each recording

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 IGM.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.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.seg - the bins in each second of IBI. Here bins.in.seg 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.

# plot burst pdfs for each recording

# write burst feature tables for each recording

printing network spikes data

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

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


# write network spike data to output file

# Check the graphs and csvs printed under the analysis$output.perDIV.dir path

Aggregating recordings

One of the strong advantages of IGM.MEA 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.

spike.features = IGM.aggregate.features(s, "spike",parameters)
ns.features = IGM.aggregate.features(s, "ns",parameters)
burst.features = IGM.aggregate.features(s, "burst",parameters)

# printing 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.

# 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.filter.maximum.DIV.inactive.ratio))

# filter network burst wells
nb.features <- lapply(nb.features, function(x) filter.wells(x, nae, parameters$well.min.rate, parameters$well.filter.maximum.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

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.

#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.

suppressMessages(permute.features.and.plot(s, "untreated", parameters$perm.n, spike.features, "spikes", analysis$output.dir))
suppressMessages(permute.features.and.plot(s, "untreated", parameters$perm.n, burst.features, "bursts", analysis$output.dir))
suppressMessages(permute.features.and.plot(s, "untreated", parameters$perm.n, ns.features, "ns", analysis$output.dir))
suppressMessages(permute.features.and.plot(s, "untreated", 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 IGM.MEA 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:

result <- suppressWarnings(dist.perm(paste0(dirname(spkListFiles[1]),"/Analysis/outputPerDIV/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_IBI_distributions.csv"),1000,"untreated","treatX"))

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:

suppressWarnings(result <- dist.perm(paste0(dirname(spkListFiles[1]),"/Analysis/outputPerDIV/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_IBI_distributions.csv"),1000,"untreated","treatX"))

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 IGM.MEA'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 IGM.MEA has to offer and use it to deeply and fully understand your MEA recordings. Please read the IGM.MEA 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! IGM MEA team

Try the IGM.MEA package in your browser

Any scripts or data that you put into this service are public.

IGM.MEA documentation built on May 29, 2017, 11:07 p.m.