Introduction to pynapple#

The goal of this tutorial is to quickly learn enough about pynapple to get started with data analysis. This tutorial assumes familiarity with the basics functionalities of numpy.

You can check how to install pynapple here.

Important

By default, pynapple will assume a time units in seconds when passing timestamps array or time parameters such as bin size (unless specified with the time_units argument)


Importing pynapple#

The convention is to import pynapple with a namespace:

import pynapple as nap

Hide code cell content

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

custom_params = {"axes.spines.right": False, "axes.spines.top": False}
sns.set_theme(style="ticks", palette="colorblind", font_scale=1.5, rc=custom_params)

Instantiating pynapple objects#

nap.Tsd: 1-dimensional time series#

If you have a 1-dimensional time series, you use the nap.Tsd object. The arguments t and d are the arguments for timestamps and data.

tsd = nap.Tsd(t=np.arange(100), d=np.random.rand(100))

print(tsd)
Time (s)
----------  ---------
0.0         0.264964
1.0         0.985347
2.0         0.0577569
3.0         0.260397
4.0         0.396245
5.0         0.487802
6.0         0.222975
...
93.0        0.931301
94.0        0.739844
95.0        0.13768
96.0        0.659098
97.0        0.0326288
98.0        0.276363
99.0        0.711838
dtype: float64, shape: (100,)

nap.TsdFrame: 2-dimensional time series#

If you have a 2-dimensional time series, you use the nap.TsdFrame object. The arguments t and d are the arguments for timestamps and data. You can add the argument columns to label each columns.

tsdframe = nap.TsdFrame(
    t=np.arange(100), d=np.random.rand(100, 3), columns=["a", "b", "c"]
)

print(tsdframe)
Time (s)           a         b          c
----------  --------  --------  ---------
0.0         0.314044  0.289517  0.579867
1.0         0.85198   0.721392  0.0992244
2.0         0.737719  0.641765  0.33217
3.0         0.606528  0.636669  0.117614
4.0         0.325479  0.704132  0.285505
5.0         0.234349  0.217875  0.649469
6.0         0.106414  0.803321  0.986147
...
93.0        0.894663  0.577591  0.763163
94.0        0.914889  0.927421  0.90118
95.0        0.472258  0.886018  0.455188
96.0        0.663691  0.456     0.446328
97.0        0.675974  0.372702  0.575793
98.0        0.150675  0.176504  0.680018
99.0        0.555286  0.153809  0.435934
dtype: float64, shape: (100, 3)

nap.TsdTensor: n-dimensional time series#

If you have larger than 2 dimensions time series (typically movies), you use the nap.TsdTensor object . The arguments t and d are the arguments for timestamps and data.

tsdtensor = nap.TsdTensor(
    t=np.arange(100), d=np.random.rand(100, 3, 4)
)

print(tsdtensor)
Time (s)
----------  -----------------------------
0.0         [[0.612788 ... 0.685036] ...]
1.0         [[0.052689 ... 0.414922] ...]
2.0         [[0.899389 ... 0.176762] ...]
3.0         [[0.298525 ... 0.634041] ...]
4.0         [[0.302505 ... 0.738074] ...]
5.0         [[0.986919 ... 0.020845] ...]
6.0         [[0.071385 ... 0.93935 ] ...]
...
93.0        [[0.159059 ... 0.464461] ...]
94.0        [[0.722966 ... 0.192899] ...]
95.0        [[0.363373 ... 0.848614] ...]
96.0        [[0.526242 ... 0.221593] ...]
97.0        [[0.166319 ... 0.828723] ...]
98.0        [[0.952259 ... 0.328062] ...]
99.0        [[0.802441 ... 0.210039] ...]
dtype: float64, shape: (100, 3, 4)

nap.IntervalSet: intervals#

The IntervalSet object stores multiple epochs with a common time unit in a table format. The epochs are strictly non-overlapping. Both start and end arguments are necessary.

epochs = nap.IntervalSet(start=[0, 10], end=[5, 15])

print(epochs)
  index    start    end
      0        0      5
      1       10     15
shape: (2, 2), time unit: sec.

nap.Ts: timestamps#

The Ts object stores timestamps data (typically spike times or reward times). The argument t for timestamps is necessary.

ts = nap.Ts(t=np.sort(np.random.uniform(0, 100, 10)))

print(ts)
Time (s)
13.205982920668157
20.772227548247958
37.9213050006308
47.34381481632922
61.65663588904818
64.50083264690902
66.249642497868
82.10355674615123
83.76865590217261
97.09773013174213
shape: 10

nap.TsGroup: group of timestamps#

TsGroup is a dictionary that stores multiple time series with different time stamps (.i.e. a group of neurons with different spike times from one session). The first argument data can be a dictionary of Ts, Tsd or numpy 1d array.

data = {
    0: nap.Ts(t=np.sort(np.random.uniform(0, 100, 1000))),
    1: nap.Ts(t=np.sort(np.random.uniform(0, 100, 2000))),
    2: nap.Ts(t=np.sort(np.random.uniform(0, 100, 3000))),
}

tsgroup = nap.TsGroup(data)

print(tsgroup, "\n")
  Index     rate
-------  -------
      0  10.0016
      1  20.0031
      2  30.0047 

Interaction between pynapple objects#

Time support : attribute of time series#

A key feature of how pynapple manipulates time series is an inherent time support object defined for Ts, Tsd, TsdFrame and TsGroup objects. The time support object is defined as an IntervalSet that provides the time series with a context. For example, the restrict operation will automatically update the time support object for the new time series. Ideally, the time support object should be defined for all time series when instantiating them. If no time series is given, the time support is inferred from the start and end of the time series.

In this example, a Tsd is instantiated with and without a time support of intervals 0 to 5 seconds. Notice how the shape of the Tsd varies.

time_support = nap.IntervalSet(start=0, end=2)

print(time_support)
  index    start    end
      0        0      2
shape: (1, 2), time unit: sec.

Without time support :

print(nap.Tsd(t=[0, 1, 2, 3, 4], d=[0, 1, 2, 3, 4]))
Time (s)
----------  --
0            0
1            1
2            2
3            3
4            4
dtype: int64, shape: (5,)

With time support :

print(
    nap.Tsd(
        t=[0, 1, 2, 3, 4], d=[0, 1, 2, 3, 4], 
        time_support = time_support
        )
    )
Time (s)
----------  --
0            0
1            1
2            2
dtype: int64, shape: (3,)

The time support object has become an attribute of the time series. Depending on the operation applied to the time series, it will be updated.

tsd = nap.Tsd(
    t=np.arange(10), d=np.random.randn(10), 
    time_support = time_support
    )

print(tsd.time_support)
  index    start    end
      0        0      2
shape: (1, 2), time unit: sec.

Restricting time series to epochs#

The central function of pynapple is the restrict method of Ts, Tsd, TsdFrame and TsGroup. The argument is an IntervalSet object. Only time points within the intervals of the IntervalSet object are returned. The time support of the time series is updated accordingly.

tsd = nap.Tsd(t=np.arange(10), d=np.random.randn(10))

ep = nap.IntervalSet(start=[0, 7], end=[3.5, 12.4])

print(ep)
  index    start    end
      0        0    3.5
      1        7   12.4
shape: (2, 2), time unit: sec.

From :

print(tsd)
Time (s)
----------  ---------
0            1.07647
1           -0.305074
2            0.706121
3           -0.283784
4            0.940177
5            0.42389
6           -0.651396
7           -0.028364
8            1.37072
9            0.577046
dtype: float64, shape: (10,)

to :

new_tsd = tsd.restrict(ep)

print(new_tsd)
Time (s)
----------  ---------
0            1.07647
1           -0.305074
2            0.706121
3           -0.283784
7           -0.028364
8            1.37072
9            0.577046
dtype: float64, shape: (7,)

Numpy & pynapple#

Pynapple relies on numpy to store the data. Pynapple objects behave very similarly to numpy and numpy functions can be applied directly

tsdtensor = nap.TsdTensor(t=np.arange(100), d=np.random.rand(100, 3, 4))

If a numpy function preserves the time axis, a pynapple object is returned.

In this example, averaging a TsdTensor along the second dimension returns a TsdFrame:

print(
    np.mean(tsdtensor, 1)
    )
Time (s)           0          1          2         3
----------  --------  ---------  ---------  --------
0.0         0.314036  0.609083   0.666829   0.133048
1.0         0.308685  0.660339   0.437743   0.180811
2.0         0.563794  0.471551   0.229721   0.489377
3.0         0.353099  0.510386   0.0827041  0.570403
4.0         0.534972  0.532929   0.567986   0.502511
5.0         0.546012  0.615079   0.618257   0.44142
6.0         0.439346  0.476122   0.588108   0.380727
...
93.0        0.355462  0.650854   0.822149   0.356729
94.0        0.43753   0.437902   0.534881   0.53155
95.0        0.429217  0.702306   0.527489   0.371155
96.0        0.300208  0.28065    0.906564   0.276186
97.0        0.402471  0.902553   0.596408   0.557511
98.0        0.385299  0.0709538  0.564957   0.423589
99.0        0.511748  0.483284   0.541785   0.55242
dtype: float64, shape: (100, 4)

Averaging along the time axis will return a numpy array object:

print(
    np.mean(tsdtensor, 0)
    )
[[0.51718591 0.53650841 0.53860603 0.48562758]
 [0.48634918 0.48249895 0.48720545 0.43578561]
 [0.49776574 0.49963704 0.48593974 0.51741589]]

Slicing objects#

Slicing time series and intervals#

Like numpy array#

Ts, Tsd, TsdFrame, TsdTensor and IntervalSet can be sliced similar to numpy array:

tsdframe = nap.TsdFrame(t=np.arange(10)/10, d=np.random.randn(10,4))
print(tsdframe)
Time (s)             0          1            2          3
----------  ----------  ---------  -----------  ---------
0            0.793046   -0.678372  -0.147684     0.420467
0.1         -0.239151   -0.644652  -1.54224      0.692377
0.2          0.388894    1.74975   -0.226464    -0.352403
0.3         -1.85676     0.174351  -0.322642     1.52344
0.4          0.282488   -0.765026  -0.587007    -0.586281
0.5          0.174353    0.509472   0.823064    -0.723782
0.6          1.191       1.03416   -1.2803      -1.80696
0.7          0.0715124  -0.924648  -0.217215     0.775025
0.8         -1.07275    -0.135049   1.89676     -1.24018
0.9         -1.40493    -1.36248   -0.00754852   1.38239
dtype: float64, shape: (10, 4)
print(tsdframe[4:7])
Time (s)           0          1          2          3
----------  --------  ---------  ---------  ---------
0.4         0.282488  -0.765026  -0.587007  -0.586281
0.5         0.174353   0.509472   0.823064  -0.723782
0.6         1.191      1.03416   -1.2803    -1.80696
dtype: float64, shape: (3, 4)
print(tsdframe[:,0])
Time (s)
----------  ----------
0            0.793046
0.1         -0.239151
0.2          0.388894
0.3         -1.85676
0.4          0.282488
0.5          0.174353
0.6          1.191
0.7          0.0715124
0.8         -1.07275
0.9         -1.40493
dtype: float64, shape: (10,)
ep = nap.IntervalSet(start=[0, 10, 20], end=[4, 15, 25])
print(ep)
  index    start    end
      0        0      4
      1       10     15
      2       20     25
shape: (3, 2), time unit: sec.
print(ep[0:2])
  index    start    end
      0        0      4
      1       10     15
shape: (2, 2), time unit: sec.
print(ep[1])
  index    start    end
      0       10     15
shape: (1, 2), time unit: sec.

Like pandas DataFrame#

Important

This page references all the way to slice TsdFrame

TsdFrame can be sliced like pandas DataFrame when the columns have been labelled with strings :

tsdframe = nap.TsdFrame(t=np.arange(10), d=np.random.randn(10,3), columns=['a', 'b', 'c'])
print(tsdframe['a'])
Time (s)
----------  ----------
0           -0.619614
1            0.334187
2            0.976863
3           -0.0250872
4           -0.815307
5            1.31271
6            0.37834
7           -1.27991
8           -0.789509
9           -0.24523
dtype: float64, shape: (10,)

but integer-indexing only works like numpy if a list of integers is used to label columns :

tsdframe = nap.TsdFrame(t=np.arange(4), d=np.random.randn(4,3), columns=[3, 2, 1])
print(tsdframe, "\n")
print(tsdframe[3])
Time (s)            3          2          1
----------  ---------  ---------  ---------
0            1.16338   -0.248389  -0.132604
1            1.59962    0.762943   0.187607
2            0.581335   0.725253   0.970056
3           -0.620604   1.22553   -1.32284
dtype: float64, shape: (4, 3) 

[-0.62060396  1.22553043 -1.32283885]

The loc method can be used to slice column-based only:

print(tsdframe.loc[3])

Slicing TsGroup#

TsGroup object can be indexed to return directly the timestamp object or sliced to return a new TsGroup.

Indexing:

print(tsgroup[0], "\n")
Time (s)
0.11775153974729058
0.16890844319858989
0.21525554185477525
0.3963934528200319
0.43807372441305725
0.4695005314430678
0.7586548750065591
...
99.70375128361518
99.71116384901578
99.75509194084225
99.76985908799779
99.78862387964787
99.80474604620979
99.93394358657973
shape: 1000 

Slicing:

print(tsgroup[[0, 2]])
  Index     rate
-------  -------
      0  10.0016
      2  30.0047

Core functions#

Objects have methods that can help transform and refine time series. This is a non exhaustive list.

Binning: counting events#

Time series objects have the count method that count the number of timestamps. This is typically used when counting number of spikes within a particular bin over multiple intervals. The returned object is a Tsd or TsdFrame with the timestamps being the center of the bins.

count = tsgroup.count(1)

print(count)
Time (s)              0    1    2
------------------  ---  ---  ---
0.5061815643453533    9   24   42
1.506181564          11   24   22
2.506181564          10   22   37
3.506181564          15   16   27
4.506181564           8   19   33
5.506181564          15   24   34
6.506181564          11    8   32
...
93.506181564         11   14   28
94.506181564         10   20   27
95.506181564         11   11   34
96.506181564         11   28   27
97.506181564         13   18   49
98.506181564         10   16   30
99.506181564         15   16   34
dtype: int64, shape: (100, 3)

Thresholding#

Some time series have specific methods. The threshold method of Tsd returns a new Tsd with all the data above or below a given value.

tsd = nap.Tsd(t=np.arange(10), d=np.random.rand(10))

print(tsd)

print(tsd.threshold(0.5))
Time (s)
----------  --------
0           0.866601
1           0.696576
2           0.624448
3           0.469884
4           0.557957
5           0.581622
6           0.966952
7           0.684131
8           0.28417
9           0.342328
dtype: float64, shape: (10,)
Time (s)
----------  --------
0           0.866601
1           0.696576
2           0.624448
4           0.557957
5           0.581622
6           0.966952
7           0.684131
dtype: float64, shape: (7,)

An important aspect of the tresholding is that the time support get updated based on the time points remaining. To get the epochs above/below a certain threshold, you can access the time support of the returned object.

print(tsd.time_support)

print(tsd.threshold(0.5, "below").time_support)
  index    start    end
      0        0      9
shape: (1, 2), time unit: sec.
  index    start    end
      0      2.5    3.5
      1      7.5    9
shape: (2, 2), time unit: sec.

Time-bin averaging of data#

Many analyses requires to bring time series to the same rates and same dimensions. A quick way to downsample a time series to match in size for example a count array is to bin average. The bin_average method takes a bin size in unit of time.

tsdframe = nap.TsdFrame(t=np.arange(0, 100)/10, d=np.random.randn(100,3))

print(tsdframe)
Time (s)             0          1           2
----------  ----------  ---------  ----------
0.0          0.407243    0.825169  -0.0132361
0.1          1.79208    -0.604754  -1.26453
0.2          0.0914083   0.575155  -0.721499
0.3         -0.2137     -1.36781   -3.29797
0.4          0.063714   -2.06433    0.0514407
0.5          0.996631   -0.82028    0.258193
0.6          0.696111    0.738647  -1.97077
...
9.3         -0.21982     0.362339   0.728657
9.4         -0.562612   -0.927121  -1.17936
9.5         -1.19439     0.263636  -1.05886
9.6          1.0972      1.92208   -0.473889
9.7          0.758236   -2.10363   -0.343735
9.8          0.156411    1.66253   -0.487307
9.9          0.217071    0.906967   1.53327
dtype: float64, shape: (100, 3)

Here we go from a timepoint every 100ms to a timepoint every second.

print(tsdframe.bin_average(1))
Time (s)             0           1           2
----------  ----------  ----------  ----------
0.5          0.460071   -0.310865   -0.784824
1.5         -0.0638389  -0.247168   -0.146775
2.5         -0.721028   -0.0347479   0.260775
3.5          0.123996    0.732123    0.288129
4.5          0.052912    0.0385481  -0.124501
5.5          0.178401    0.170248   -0.222094
6.5          0.466871    0.0560009   0.0381009
7.5          0.311081   -0.0525766   0.0490344
8.5          0.131658    0.639697   -0.211918
9.5         -0.162995   -0.0930361  -0.382216
dtype: float64, shape: (10, 3)

Loading data#

See here for more details about loading data.

Loading NWB#

Pynapple supports by default the NWB standard.

NWB files can be loaded with :

nwb = nap.load_file("path/to/my.nwb")

or directly with the NWBFile class:

nwb = nap.NWBFile("path/to/my.nwb")

print(nwb)
my.nwb
┍━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━┑
│ Keys            │ Type        │
┝━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━┥
│ units           │ TsGroup     │
┕━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━┙

The returned object behaves like a dictionary. The first column indicates the keys. The second column indicate the object type.

print(nwb['units'])
  Index    rate  location      group
-------  ------  ----------  -------
      0    1.0  brain        0
      1    1.0  brain        0
      2    1.0  brain        0

Overview of advanced analysis#

The process module of pynapple contains submodules that group methods that can be applied for high level analysis. All of the method are directly available from the nap namespace.

Discrete correlograms & ISI

This module analyses discrete events, specifically correlograms (for example by computing the cross-correlograms of a population of neurons) and interspike interval (ISI) distributions.

Bayesian decoding

The decoding module performs bayesian decoding given a set of tuning curves and a TsGroup.

Filtering

Bandpass, lowpass, highpass or bandstop filtering can be done to any time series using either Butterworth filter or windowed-sinc convolution.

Perievent time histogram

The perievent module has a set of functions to center time series and timestamps data around a particular events.

Randomizing

The randomize module holds multiple technique to shuffle timestamps in order to create surrogate datasets.

Spectrum

The spectrum module contains the methods to return the (mean) power spectral density of a time series.

Tuning curves

Tuning curves of neurons based on spiking or calcium activity can be computed.

Wavelets

The wavelets module performs Morlet wavelets decomposition of a time series.

Phases & envelopes

This modules allows for computing analytic signals and extracting the phase and envelope.

Warping

This module provides methods for building trial-based tensors and time-warped trial-based tensors.