Plot spikes and network bursts#

Generate images from datasets and files from previous analyses containing spikes and network bursts.#

autoMEA can also be used to visualize activity from a dataset that was previously analyzed - for which spikes and bursts were already detected.

Here we will show how to, given a dataset containing signal, spikes and bursts, how to generate images to visualize the analyzed data, and compare the results with those obtained using autoMEA.

We first import autoMEA, together with numpy and pandas, which will be used to manipulate data.

import automea
import numpy as np 
import pandas as pd 
am = automea.Analysis()

We create a variable containing the path to the dataset, the dataset name, and the well we want to analyze in the dataset.

path_to_dataset = '/Volumes/T7/automea_files_to_plot/'
dataset = 'MJ22C002-1782_DIV7.h5'
well = 'D5'

And we create separate folders to save signals and spikes rasters plots.

import os  
plots_dir = f'plots/{dataset}_{well}'
if os.path.exists(plots_dir) is False:
    os.mkdir(plots_dir)
    os.mkdir(plots_dir+'/signals')
    os.mkdir(plots_dir+'/rasters')

We inform autoMEA the path to find the dataset, and load it using the loadh5() method.

am.path_to_dataset = path_to_dataset
am.loadh5(dataset)

After loaded, we can check that autoMEA is already storing a signal array.

am.signal, am.signal.shape

(array([[ 27,  23, -17, ..., -24,  31,  31],
        [-70,  15,  69, ...,  42,   4,  -8],
        [ -5, -25, -41, ..., -13, -54, -11],
        ...,
        [-17, -14,   5, ..., -17, -80,  17],
        [-36,   7,  -7, ..., -37, -19, -29],
        [  8, -20,  22, ...,   9,  73,  71]], dtype=int32),
 (237, 6000000))

Following, we use the loadsignal() method, to keep in memory only the signal from the well we want to analyze.

am.loadsignal(am.signal[am.wellsFromData == am.wellLabelIndexDict[well]])

Here we load a model, and use autoMEA to detect all the variables we want to compare with the already-analyzed dataset.

am.model_name = 'signal30.h5'
am.loadmodel()

Metal device set to: Apple M1

systemMemory: 8.00 GB
maxCacheSize: 2.67 GB

am.detect_threshold()
am.detect_spikes()
am.detect_reverbs(method = 'model')
am.detect_bursts()
am.detect_net('reverbs')
am.detect_net('bursts')

Now we load the previously-analyzed spikes and network bursts, from deticated ‘csv’ files. We filter the dataset to keep only the activity from the well we want to analyze, and save the spikes in an array named spikes_manual, and the network bursts in a net_manual array.

path_to_manual_files = '/Volumes/T7/automea_files_to_plot/Examples list/RHEB_DIV07/'
spikes_manual_file = 'Timestamps_MJ22C002_1782_DIV07_D5.csv'
netbursts_manual_file = 'NetworkBursts_MJ22C002_1782_DIV07_D5.csv'

spikes_manual_df = pd.read_csv(path_to_manual_files + spikes_manual_file, encoding='unicode_escape')

spikes_manual_df

	Channel ID	Channel Label	Well ID	Well Label	Compound ID	Compound Name	Experiment	Dose [pM]	Dose Label	Timestamp [µs]	Maximum Amplitude [pV]	Minimum Amplitude [pV]	Peak-to-peak Amplitude [pV]
0	191	21	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	426300	16212288	-25391304	41603592
1	191	21	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	435100	5066340	-11920800	16987140
2	191	21	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	452400	7569708	-14006940	21576648
3	191	21	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	494500	7808124	-12516840	20324964
4	191	21	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	503100	50007756	-59246376	109254132
...	...	...	...	...	...	...	...	...	...	...	...	...	...
65010	185	34	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	598999200	7629312	-23185956	30815268
65011	185	34	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	599003300	15079812	-2264952	17344764
65012	185	34	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	599010600	13470504	-8702184	22172688
65013	185	34	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	599039800	16391100	-7092876	23483976
65014	185	34	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	599148400	5721984	-13887732	19609716

65015 rows × 13 columns

channel_labels = np.unique(spikes_manual_df[spikes_manual_df['Well Label'] == well]['Channel Label'].iloc[:].values)
channel_labels

array([12, 13, 21, 22, 23, 24, 31, 32, 33, 34, 42, 43])

spikes_manual = []
for channel_label in channel_labels:
    spikes_manual.append([])
    for index in spikes_manual_df[(spikes_manual_df['Channel Label'] == channel_label) & (spikes_manual_df['Well Label'] == well)].index:
        spikes_manual[-1].append(spikes_manual_df['Timestamp [µs]'][index]//100)

netbursts_manual_df = pd.read_csv(path_to_manual_files + netbursts_manual_file, encoding='unicode_escape')
netbursts_manual_df

	Well ID	Well Label	Compound ID	Compound Name	Experiment	Dose [pM]	Dose Label	Start timestamp [µs]	Duration [µs]	Spike Count	Spike Frequency [Hz]
0	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	494500	250600	296	98.430
1	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	14255400	319800	608	158.432
2	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	18646700	444100	749	140.546
3	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	26212300	305800	653	177.949
4	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	30838400	265700	603	189.123
...	...	...	...	...	...	...	...	...	...	...	...
91	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	573895000	311100	216	57.859
92	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	578123900	375300	725	160.982
93	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	585406500	363400	487	111.677
94	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	591906400	390700	690	147.172
95	22	D5	No compound	No compound	20181127_No_compound_m201118_DIV7_rescue4EB_001	0	Control	598675200	379000	734	161.390

96 rows × 11 columns

net_manual = []
for index in netbursts_manual_df[netbursts_manual_df['Well Label']==well].index:
    start_ts = netbursts_manual_df['Start timestamp [µs]'][index]
    duration = netbursts_manual_df['Duration [µs]'][index]
    net_manual.append([start_ts//100, (start_ts+duration)//100])

After having spikes and network bursts both from a previously performed analysis, and obtained with autoMEA, we use the plot_raster() method to visualize raster plots, so that the two version can be compared.

duration = 60
duration_shorter = 5

for i, start_time in enumerate(np.random.randint(0,500,size=3)):

    end_time = start_time + duration

    am.plot_raster(am.spikes, net_bursts = am.net_reverbs, start_time = start_time, end_time = end_time, 
                   save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_allwell_prednetreverbs.eps', show = False)
    am.plot_raster(am.spikes, net_bursts = net_manual, start_time = start_time, end_time = end_time, 
                   save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_allwell_manual.eps', show = False)

    for j in range(5):
        start_time_shorter = start_time + np.random.randint(0,95)
        channel_random  = np.random.randint(len(am.signal))
        end_time_shorter = start_time_shorter + duration_shorter
        am.plot_raster(am.spikes[channel_random], net_bursts = am.net_reverbs, start_time = start_time_shorter, end_time = end_time_shorter, 
                        save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_zoomedstart{start_time+start_time_shorter}_zoomedduration{duration_shorter}_prednetreverbs.eps', show = False)
        am.plot_raster(am.spikes[channel_random], net_bursts = net_manual, start_time = start_time_shorter, end_time = end_time_shorter,
                        save = f'{plots_dir}/rasters/raster_start{start_time}_duration{duration}_zoomedstart{start_time+start_time_shorter}_zoomedduration{duration_shorter}_manual.eps', show = False)

And finally we plot a (randomly-selected) 5-second activity with signal and threshold, together with autoMEA and manual network bursts, to compare them.

duration_shorter = 5

for _ in range(20):
    random_net = am.net_reverbs[np.random.randint(len(am.net_reverbs))]
    middle_of_net = random_net[0] + ((random_net[-1] - random_net[0])/2)
    start_time = middle_of_net/am.samplingFreq - 2.5
    channel_random  = np.random.randint(len(am.signal))


    am.plot_window(am.signal[channel_random], start_time=start_time, duration=duration_shorter,
                    threshold=am.threshold[channel_random], spikes = am.spikes[channel_random],
                    bursts=am.net_reverbs, net_bursts=net_manual, yunits = 'V', show = False,
                    save = f'{plots_dir}/signals/signal_channel{channel_random}_start{start_time}_duration5.eps')

Plot spikes and network bursts

Contents

Plot spikes and network bursts#

Generate images from datasets and files from previous analyses containing spikes and network bursts.#