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.
The IGM.MEA package is available for installation under the CRAN repository:
#install.packages( "IGM.MEA",repos="http://cran.us.r-project.org")
To load the package please load IGM.MEA and dependencies:
library(IGM.MEA)
## Warning: no DISPLAY variable so Tk is not available
library(plyr)
library(ggplot2)
library(reshape2)
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")
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
data.dir<-dirname(spkListFiles[1])
# 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)
}
## Total number of spikes: 102439
## Unique number of electrodes: 255
## Time range [0.010 35.016] (seconds)
## nelectrodes nspikes time.min
## exampleRecording_1012016_plate1_DIV1.gz 255 102439 0.01
## time.max
## exampleRecording_1012016_plate1_DIV1.gz 35.01592
## [1] "spikes" "sCount" "epos" "names" "array"
## [6] "treatment" "dose" "size" "well" "units"
## Total number of spikes: 84758
## Unique number of electrodes: 268
## Time range [0.010 50.306] (seconds)
## nelectrodes nspikes time.min
## exampleRecording_1012016_plate1_DIV3.gz 268 84758 0.01048
## time.max
## exampleRecording_1012016_plate1_DIV3.gz 50.30616
## [1] "spikes" "sCount" "epos" "names" "array"
## [6] "treatment" "dose" "size" "well" "units"
## Total number of spikes: 84021
## Unique number of electrodes: 252
## Time range [0.010 50.058] (seconds)
## nelectrodes nspikes time.min
## exampleRecording_1012016_plate1_DIV4.gz 252 84021 0.01
## time.max
## exampleRecording_1012016_plate1_DIV4.gz 50.05832
## [1] "spikes" "sCount" "epos" "names" "array"
## [6] "treatment" "dose" "size" "well" "units"
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.
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
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
parameters$burst.type="ps"
# Construct the 'spike.list' object and calculate spike features
s<-calculate.spike.features(analysis$Robject, parameters)
# 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 <- 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:
s[[1]]$spikes$B3_41
## [1] 0.55344 2.84256 9.51040 10.51240 12.84232 14.17480 20.87056
## [8] 21.53936 23.87976 26.23032 30.24168
To view burst information for bursts calculated for electrode E7_42 of the 2nd recording, try the following:
s[[2]]$allb$E7_42
## beg end IBI len durn mean.isis SI
## [1,] 11 16 NA 6 0.55776 0.1115520 5.532836
## [2,] 29 31 8.79088 3 0.03752 0.0187600 6.099738
## [3,] 39 43 5.86656 5 0.31488 0.0787200 5.846385
## [4,] 57 63 14.88784 7 0.60760 0.1012667 6.895707
## [5,] 64 80 5.17472 17 2.63536 0.1647100 10.018063
## [6,] 84 86 3.79088 3 0.05944 0.0297200 5.206044
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.
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
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 :
s[[i]]$ns.all$B5$en.brief
## spikes ns spikes.in.ns percent.of.spikes.in.ns mean.spikes.per.ns
## B5_11 86 2 2 2.325581 1.000000
## B5_13 103 2 2 1.941748 1.000000
## B5_21 22 3 4 18.181818 1.333333
## B5_23 107 4 4 3.738318 1.000000
## B5_31 17 4 4 23.529412 1.000000
## B5_32 5 0 0 0.000000 NA
## B5_34 104 4 4 3.846154 1.000000
## B5_41 11 2 2 18.181818 1.000000
## B5_42 19 1 1 5.263158 1.000000
## B5_44 47 3 4 8.510638 1.333333
## sd.spikes.per.ns mean.insis
## B5_11 0.0000000 7.700000
## B5_13 0.0000000 0.720000
## B5_21 0.5773503 3.805000
## B5_23 0.0000000 2.566667
## B5_31 0.0000000 2.536667
## B5_32 NA NA
## B5_34 0.0000000 1.096667
## B5_41 0.0000000 4.410000
## B5_42 NA NA
## B5_44 0.5773503 1.645000
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,
parameters$min_electrodes,
parameters$local_region_min_nAE)
## calculating network bursts for recording exampleRecording_1012016_plate1_DIV1.RData
## calculating network bursts for recording exampleRecording_1012016_plate1_DIV3.RData
## Warning in min(well.data[well.data > 0]): no non-missing arguments to min;
## returning Inf
## calculating network bursts for recording exampleRecording_1012016_plate1_DIV4.RData
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]]
}
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.
# print spikes graphs (pdf format) for each recording
IGM.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))
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
suppressWarnings(IGM.plot.plate.summary.for.bursts(s,analysis$output.perDIV.dir,parameters))
## [1] "Running IBI distribution analysis."
## Arguments: minVals=15; xlimit=20; binsInSec=5;perWell=0; duration=35.00592; feature='IBI'; filterValuesByMin=0; minValues=0
## [1] "Running ISI distribution analysis."
## Arguments: minVals=15; xlimit=0.5; binsInSec=100;perWell=0; duration=35.00592; feature='ISI'; filterValuesByMin=0; minValues=0
## [1] "Running nSpikes in bursts distribution analysis."
## Arguments: minVals=5; xlimit=200; binsInSec=1;perWell=0; duration=35.00592; feature='nspikesInBurst'; filterValuesByMin=0; minValues=0
## [1] "Running duration of bursts distribution analysis."
## Arguments: minVals=15; xlimit=18; binsInSec=10;perWell=0; duration=35.00592; feature='duration'; filterValuesByMin=0; minValues=0
## [1] "Running spike density in bursts distribution analysis."
## Arguments: minVals=15; xlimit=300; binsInSec=1;perWell=0; duration=35.00592; feature='spikesDensityInBurst'; filterValuesByMin=0; minValues=0
## [1] "Running IBI distribution analysis."
## Arguments: minVals=15; xlimit=20; binsInSec=5;perWell=0; duration=50.29568; feature='IBI'; filterValuesByMin=0; minValues=0
## [1] "Running ISI distribution analysis."
## Arguments: minVals=15; xlimit=0.5; binsInSec=100;perWell=0; duration=50.29568; feature='ISI'; filterValuesByMin=0; minValues=0
## [1] "Running nSpikes in bursts distribution analysis."
## Arguments: minVals=5; xlimit=200; binsInSec=1;perWell=0; duration=50.29568; feature='nspikesInBurst'; filterValuesByMin=0; minValues=0
## [1] "Running duration of bursts distribution analysis."
## Arguments: minVals=15; xlimit=18; binsInSec=10;perWell=0; duration=50.29568; feature='duration'; filterValuesByMin=0; minValues=0
## [1] "Running spike density in bursts distribution analysis."
## Arguments: minVals=15; xlimit=300; binsInSec=1;perWell=0; duration=50.29568; feature='spikesDensityInBurst'; filterValuesByMin=0; minValues=0
## [1] "Running IBI distribution analysis."
## Arguments: minVals=15; xlimit=20; binsInSec=5;perWell=0; duration=50.04832; feature='IBI'; filterValuesByMin=0; minValues=0
## [1] "Running ISI distribution analysis."
## Arguments: minVals=15; xlimit=0.5; binsInSec=100;perWell=0; duration=50.04832; feature='ISI'; filterValuesByMin=0; minValues=0
## [1] "Running nSpikes in bursts distribution analysis."
## Arguments: minVals=5; xlimit=200; binsInSec=1;perWell=0; duration=50.04832; feature='nspikesInBurst'; filterValuesByMin=0; minValues=0
## [1] "Running duration of bursts distribution analysis."
## Arguments: minVals=15; xlimit=18; binsInSec=10;perWell=0; duration=50.04832; feature='duration'; filterValuesByMin=0; minValues=0
## [1] "Running spike density in bursts distribution analysis."
## Arguments: minVals=15; xlimit=300; binsInSec=1;perWell=0; duration=50.04832; feature='spikesDensityInBurst'; filterValuesByMin=0; minValues=0
# write burst feature tables for each recording
write.plate.summary.for.bursts(s,analysis$output.perDIV.dir)
The commands below print single recording ns data. Here they are printed for the first of the three recordings
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"))
IGM.xyplot.network.spikes(nspikes)
IGM.plot.active.wells.network.spikes(nspikes)
dev.off()
## pdf
## 2
# 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
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
spike.features$nAE
## well treatment DIV1 DIV3 DIV4
## 1 B2 treatX 14 13 13
## 2 B3 treatX 16 16 16
## 3 B4 untreated 13 13 <NA>
## 4 B5 untreated 10 <NA> 10
## 5 B6 treatY 11 12 12
## 6 B7 treatY 12 12 12
## 7 C2 treatX 14 13 14
## 8 C3 treatX 13 13 13
## 9 C6 treatY 10 12 10
## 10 C7 treatY 12 13 13
## 11 D2 treatX 15 14 14
## 12 D3 treatX 16 15 15
## 13 D5 untreated <NA> 14 <NA>
## 14 D6 treatY 11 13 13
## 15 D7 treatY 13 12 12
## 16 E2 treatX 16 14 15
## 17 E3 treatX 16 16 15
## 18 E4 untreated <NA> 14 15
## 19 E5 untreated 12 11 12
## 20 E6 treatY 14 14 13
## 21 E7 treatY 14 14 14
#Feel free to explore the spike/ns/burst and nb.features for the different features they offer
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
spike.features$nAE
## well treatment DIV1 DIV3 DIV4
## 1 B2 treatX 14 13 13
## 2 B3 treatX 16 16 16
## 3 B4 untreated 13 13 0
## 4 B5 untreated 10 0 10
## 5 B6 treatY 11 12 12
## 6 B7 treatY 12 12 12
## 7 C2 treatX 14 13 14
## 8 C3 treatX 13 13 13
## 9 C6 treatY 10 12 10
## 10 C7 treatY 12 13 13
## 11 D2 treatX 15 14 14
## 12 D3 treatX 16 15 15
## 14 D6 treatY 11 13 13
## 15 D7 treatY 13 12 12
## 16 E2 treatX 16 14 15
## 17 E3 treatX 16 16 15
## 18 E4 untreated 0 14 15
## 19 E5 untreated 12 11 12
## 20 E6 treatY 14 14 13
## 21 E7 treatY 14 14 14
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).
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")
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:
result <- suppressWarnings(dist.perm(paste0(dirname(spkListFiles[1]),"/Analysis/outputPerDIV/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_IBI_distributions.csv"),1000,"untreated","treatX"))
## 100 permutations
## 200 permutations
## 300 permutations
## 400 permutations
## 500 permutations
## 600 permutations
## 700 permutations
## 800 permutations
## 900 permutations
## 1000 permutations
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:
suppressWarnings(result <- dist.perm(paste0(dirname(spkListFiles[1]),"/Analysis/outputPerDIV/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_IBI_distributions.csv"),1000,"untreated","treatX"))
## 100 permutations
## 200 permutations
## 300 permutations
## 400 permutations
## 500 permutations
## 600 permutations
## 700 permutations
## 800 permutations
## 900 permutations
## 1000 permutations
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 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