# meaRtools: Performing a Comprehensive Analysis of Multiple MEA Recordings

#### Sahar Gelfman

#### 2018-07-06

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 two sequential recordings of the same plate.

## Installing the package

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

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

## Loading the package

To load the package please load meaRtools and dependencies:

    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/ :

    # set path to "_spike_list.csv" files from the file path in 'filesPath'
    spk_list_files<-c(
                    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

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

    ## [1] "Creating output folders under /var/folders/mn/tjk6zsyd6cz3rxb5vwkyn0p80000gp/T//RtmpgZhGda"

    # 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

    # A loop to go over all 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)
    }

    ## Total number of spikes: 100957
    ## Unique number of electrodes: 269
    ## Time range [0.010 57.678] (seconds)
    ##                                         nelectrodes nspikes time.min
    ## exampleRecording_1012016_plate1_DIV3.gz         269  100957  0.01048
    ##                                         time.max
    ## exampleRecording_1012016_plate1_DIV3.gz 57.67824
    ##  [1] "spikes"    "scount"    "epos"      "names"     "array"
    ##  [6] "treatment" "dose"      "size"      "well"      "units"
    ## Total number of spikes: 106530
    ## Unique number of electrodes: 268
    ## Time range [0.010 61.761] (seconds)
    ##                                         nelectrodes nspikes time.min
    ## exampleRecording_1012016_plate1_DIV4.gz         268  106530     0.01
    ##                                         time.max
    ## exampleRecording_1012016_plate1_DIV4.gz  61.7608
    ##  [1] "spikes"    "scount"    "epos"      "names"     "array"
    ##  [6] "treatment" "dose"      "size"      "well"      "units"

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

    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.

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

    s[[1]]$spikes$B3_41

    ## [1] 14.83456

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

    s[[1]]$allb$E7_42

    ##      beg end      ibi len    durn mean_isis si
    ## [1,]  11  16       NA   6 0.55776 0.1115520  1
    ## [2,]  39  43 14.34256   5 0.31488 0.0787200  1
    ## [3,]  57  63 14.28024   7 0.60760 0.1012667  1

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

#### 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,
                                           parameters$min_electrodes,
                                           parameters$local_region_min_nae)

    ## calculating network bursts for recording exampleRecording_1012016_plate1_DIV3.RData
    ## 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]]
        }

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

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

    unlist(s[[1]]$mean_sttc)

    ##        E2        E3        E4        E5        E6        E7        D2
    ## 0.6872915 0.6679063 0.7111581 0.7648321 0.4163986 0.1673084 0.7511649
    ##        D3        D5        D6        D7        C2        C3        C6
    ## 0.6838312 0.6438838 0.2251951 0.2401321 0.7418769 0.7098640 0.3933835
    ##        C7        B2        B3        B4        B6        B7
    ## 0.1724659 0.7611882 0.7397869 0.6627209 0.4013349 0.1410473

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

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

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

    ## [1] "Running IBI distribution analysis."
    ## Arguments: min_vals=15; xlimit=20; bins_in_sec=5;per_well=0; duration=57.66776; feature='ibi'; filter_values_by_min=0; min_values=0
    ## [1] "Running ISI distribution analysis."
    ## Arguments: min_vals=15; xlimit=0.5; bins_in_sec=100;per_well=0; duration=57.66776; feature='isi'; filter_values_by_min=0; min_values=0

    ## [1] "Running nSpikes in bursts distribution analysis."
    ## Arguments: min_vals=5; xlimit=200; bins_in_sec=1;per_well=0; duration=57.66776; feature='nspikes_in_burst'; filter_values_by_min=0; min_values=0

    ## [1] "Running duration of bursts distribution analysis."
    ## Arguments: min_vals=15; xlimit=18; bins_in_sec=5;per_well=0; duration=57.66776; feature='duration'; filter_values_by_min=0; min_values=0
    ## [1] "Running spike density in bursts distribution analysis."
    ## Arguments: min_vals=15; xlimit=300; bins_in_sec=1;per_well=0; duration=57.66776; feature='spikes_density_in_burst'; filter_values_by_min=0; min_values=0

    ## [1] "Running IBI distribution analysis."
    ## Arguments: min_vals=15; xlimit=20; bins_in_sec=5;per_well=0; duration=61.7508; feature='ibi'; filter_values_by_min=0; min_values=0
    ## [1] "Running ISI distribution analysis."
    ## Arguments: min_vals=15; xlimit=0.5; bins_in_sec=100;per_well=0; duration=61.7508; feature='isi'; filter_values_by_min=0; min_values=0

    ## [1] "Running nSpikes in bursts distribution analysis."
    ## Arguments: min_vals=5; xlimit=200; bins_in_sec=1;per_well=0; duration=61.7508; feature='nspikes_in_burst'; filter_values_by_min=0; min_values=0

    ## [1] "Running duration of bursts distribution analysis."
    ## Arguments: min_vals=15; xlimit=18; bins_in_sec=5;per_well=0; duration=61.7508; feature='duration'; filter_values_by_min=0; min_values=0
    ## [1] "Running spike density in bursts distribution analysis."
    ## Arguments: min_vals=15; xlimit=300; bins_in_sec=1;per_well=0; duration=61.7508; feature='spikes_density_in_burst'; filter_values_by_min=0; min_values=0

    # 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

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

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

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

    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

    ##    well treatment div3 div4
    ## 1    B2    treatX   13   13
    ## 2    B3    treatX   16   15
    ## 3    B4 untreated   13   12
    ## 4    B5 untreated <NA>   10
    ## 5    B6    treatY   12   12
    ## 6    B7    treatY   12   12
    ## 7    C2    treatX   13   14
    ## 8    C3    treatX   13   13
    ## 9    C6    treatY   12   10
    ## 10   C7    treatY   13   13
    ## 11   D2    treatX   14   13
    ## 12   D3    treatX   15   15
    ## 13   D5 untreated   14 <NA>
    ## 14   D6    treatY   13   12
    ## 15   D7    treatY   12   12
    ## 16   E2    treatX   15   14
    ## 17   E3    treatX   16   15
    ## 18   E4 untreated   14   15
    ## 19   E5 untreated   11   12
    ## 20   E6    treatY   14   13
    ## 21   E7    treatY   14   14

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

    ##    well treatment div3 div4
    ## 1    B2    treatX   13   13
    ## 2    B3    treatX   16   15
    ## 3    B4 untreated   13   12
    ## 5    B6    treatY   12   12
    ## 6    B7    treatY   12   12
    ## 7    C2    treatX   13   14
    ## 8    C3    treatX   13   13
    ## 9    C6    treatY   12   10
    ## 10   C7    treatY   13   13
    ## 11   D2    treatX   14   13
    ## 12   D3    treatX   15   15
    ## 14   D6    treatY   13   12
    ## 15   D7    treatY   12   12
    ## 16   E2    treatX   15   14
    ## 17   E3    treatX   16   15
    ## 18   E4 untreated   14   15
    ## 19   E5 untreated   11   12
    ## 20   E6    treatY   14   13
    ## 21   E7    treatY   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).

## 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, "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 two 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:

    result <- suppressWarnings(dist_perm(paste0(analysis$output_perDIV_dir,"/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_isi_distributions.csv"),1000,"treatX","treatY"))

    ## 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="ISI")
    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(analysis$output_perDIV_dir,"/distributionFiles/exampleRecording_1012016_plate1_DATE_TIME_isi_distributions.csv"),1000,"treatX","treatY"))

    ## 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="ISI")
    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 two sequential 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 two 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