### Sidebar

Reference

analysis:nsb2016:week3long

This is an old revision of the document!

## Module 3: Visualizing neural data in MATLAB

Goals:

• Master the basics of plotting data in MATLAB using the versatile plot() function
• Understand how to use handles to fine-tune formatting
• Meet some existing visualization tools: MultiRaster(), PlotTSDfromIV() and ft_databrowser()
• Use interactive figure properties with callback functions for dynamic data viewing
• Learn to use functions with variable numbers of input arguments (using varargin and cfg inputs) effectively
• Understand which file format to export figures to
• Advanced bonus section: GUIs, movies, and sonification of spiking data

Resources:

## Introductory remarks

Visualization is an essential component of any data analysis workflow. You would never perform brain surgery or another intricate experimental procedure without a good view of what you are doing – if conditions are met for a next step, whether a manipulation has the desired outcome – and likewise, you should not do data analysis blind. This means that after every step you want to be able to see what changed.

Also, for many analyses, a visualization is part of the final outcome. Most neuroscience papers have complex multipanel figures, such as this one from Dragoi and Tonegawa (Nature, 2011):

Notice that panel Aa, a highly selective and organized view of raw data, already contains LFP data (top) and spike data from 18 cells!

This sort of figure is constructed by first producing the individual panels in MATLAB, and then perhaps combining them in a graphic design package such as Adobe Illustrator/InDesign. We will do the MATLAB part here and aim to produce something similar to panel Aa, while introducing the fundamentals of plotting in MATLAB. Other types of plots, such as bar plots, errorbars, color bitmaps, filled shapes, and shading will be introduced as they are needed in later modules.

An important principle in graphical design for scientific communication is to only include those elements which contribute to the messages(s) you want to communicate. Remove and simplify as much as possible: this will make your figures more readable. The Tufte book linked to above is a beautiful statement of this (and other) guiding principles of graphic design.

Even if you do not go as far as reading it, do at least make the labels in your figures big enough to be readable!

If you haven't done so already, this would be a good time for a git pull to make sure you have the latest version of the course codebase!

Now, make sure you have the folder R042-2013-08-18 from the data share, and that this is placed in a sensible location (NOT in a GitHub or project folder!).

This data set was collected from the dorsal CA1 region of the hippocampus of a rat performing a T-maze task, using a drive with 16 independently movable tetrodes (by Alyssa Carey for her Master's thesis work). Spike and LFP data was recorded from each tetrode; possible spike events were detected online and stored for offline spike-sorting, and LFPs were sampled at 2kHz and bandpass filtered between 1-475Hz. (A quirk of this particular data set is that certain time intervals are cut out of the spike data, but not the LFP. So you may notice some odd looking gaps in the rasterplot later.)

As in previous modules, use the shortcut button created in the previous module to set your path. In your project folder, create a folder with today's date and cd to it. This will be your working directory with stuff that you don't necessarily intend to push to GitHub.

Create a sandbox.m file in your daily folder. In this file, use cell mode to load some spike trains, a LFP, and position data as follows (recall you can use Ctrl+Enter to execute the code in a cell):

%% cd to data folder -- replace this with yours
fd = 'D:\data\DataAnalysisTutorial\R042-2013-08-18';
cd(fd);

%% load the data (note, may need to unzip position data first)
cfg = [];
S = LoadSpikes(cfg)

This should give a data structure containing spike train data:

LoadSpikes: Loading 67 files...

S =

type: 'ts'
t: {1x67 cell}
label: {1x67 cell}
cfg: [1x1 struct]
usr: [1x1 struct]

The details of this ts structure are discussed in Module 2. In brief, each cell of S.t contains the spike times from a putative neuron. The qualifier “putative” is used because this is extracellular data and spike-sorting is not perfect, so it's likely there will be some spikes missing and some spikes included that are not from this neuron. Always remember this even if I will omit the “putative” from now on for short. You can see we have loaded 67 neurons for now.

The LoadSpikes() function has some features not discussed here, such as the contents of the usr field; you can type help LoadSpikes() to learn more.

cfg = [];
cfg.fc = {'R042-2013-08-18-CSC03a.ncs'};
csc = LoadCSC(cfg)

This gives:

LoadCSC: Loading 1 file(s)...

csc =

type: 'tsd'
tvec: [8802816x1 double]
data: [1x8802816 double]
label: {'R042-2013-08-18-CSC03a.ncs'}
cfg: [1x1 struct]

csc is the common Neuralynx designation for “continuously sampled channel” and typically is an EEG or LFP type signal sampled and filtered so that high-frequency components such as spikes are not accessible. It is possible to have wide-band, 32kHz CSCs suitable for spike extraction, but these are not included in the current dataset. As discussed more extensively in Module 2 a LFP is defined by matching sets of sample timestamps (csc.tvec) and sampled data (csc.data).

Before we proceed, let's restrict the data we have loaded (about 90 minutes' worth) to a more manageable size.

A useful function that works on ts, tsd and iv objects is restrict().

☛ Use MATLAB's help function to look up the usage information for restrict(). In a new editor cell, use this function to create S_r and csc_r variables, containing the spike and LFP data respectively, restricted between 5900 and 6000 seconds.

Hint: make a habit of using variables instead of hard-coding specific values. Applying this good programming practice to the above means that you first define the interval of interest, and then use the same resulting variable to two calls to restrict(). That way, if you ever want to change this interval, you can do so in a single, clearly visible location at the beginning of your code. This is much more robust than having to scan through the whole script to find and modify each instance of some hard-coded numbers!

☛ Verify that the restriction worked by inspecting the csc_r variable – you should notice that the length of the tvec and data fields has been suitably reduced.

## Do-it-yourself plotting

If you are already familiar with MATLAB's basic plotting functions and how to use object handles to set properties of drawing objects, you can skip this section.

### Basic plot commands

Let's look at the spike times of the 17th cell in our data set (note the use of the curly brackets {}, if this is confusing, review the MATLAB help on cell arrays):

%%
iC= 17;
this_spk = S_r.t{iC}

this_spk =

1.0e+03 *

5.9046
5.9046
5.9046
5.9046

That's not a very helpful format, it looks like all four spikes have the same spike time, but note the 1.0e+03 (meaning 10^3, or 1000)!

☛ Type format bank and try again.

You should now see somewhat more clearly the spike times of the five spikes this neuron emitted (in our restriction interval).

An alternative to format is to use the ubiquitous fprintf() command:

>> fprintf('%4.3f\n',this_spk)
5904.630
5904.635
5904.638
5904.643

This way, you have precise control over the formatting of command line output – the cryptic %4.3f tag specifies that the contents of this_spk should be formatted as a floating-point number with 4 digits before, and 3 digits following, the decimal point; \n specifies a newline. (For reference: fprintf() has many other formatting options).

Now that we know what we are dealing with, let's meet the plot() command:

plot(this_spk)

You should get:

This may not be what you were expecting! If you give plot() only one input argument (this_spk in this case) it will, by default, plot the index of each element against its value.

Let's unpack that statement so it is really clear: we have an array with five values here: the spike times of five spikes of neuron 11. plot() plots the value of the first element (this_spk(1)) at x-coordinate 1 (its index, i.e. position in the array), and so forth for the whole length of the array. It also connects the data points with a blue line.

Let's plot the same data differently:

plot(this_spk,0,'.k');

This gives:

A pretty different result! Because we have given plot() multiple input arguments, it interprets the first (this_spk) as the x-coordinate, and the second (0) as the y-coordinate. Indeed you can see in the plot that the spike times now show up on the x-axis, in contrast to the first plot we made, and the y-coordinate is 0 for all spikes.

The final argument ('.k') specifies we want the data plotted as dots (.) and in black ('k').

Plotting spikes as dots is okay, but tickmarks (vertical lines, as in the example figure at the top) are better. Here is how to do it:

h = plot([this_spk this_spk],[0.5 1.5],'k');
axis([xlim 0 5])

To understand what happened here, recall that plot() uses the first argument as x-coordinates, and the second as y-coordinates. So the above plot command used our spike times (this_spk) as the x-coordinates, drawing a black line at each spike time between the specified y-coordinates [0.5 1.5].

The axis command sets the axis limit according to the 4-dimensional array [xmin xmax ymin ymax]; xlim returns the current xmin and xmax. What the h output argument is for will be explained in the next section.

☛ (Test your MATLAB skills) Extend the example code above to make a rasterplot containing the spikes for all cells in S_r. Each cell should have its spikes appear on the corresponding row, as in the example plot at the top. To facilitate re-use later, wrap your code in a function, PlotSpikeRaster, that plots the raster for any ts input. Thus, it should also work for task events, for instance, by virtue of following the common ts datatype specification. Hint: how will you handle the case where some cell doesn't have any spikes?

### Using handles to fine-tune plotting

At the end of the previous section, I specified an output argument h with the plot() command. This output is a handle, a MATLAB data type that contains various properties associated with an object. In this case, h is a handle for the spikes we just plotted, and we can use it to change their appearance.

☛ Below are a few examples of how to use a plot handle. Try them – one at a time – and see what happens!

set(h,'Color',[0 1 0]); % "set Color property of handle h to value [0 1 0] i.e. green in [red green blue] format
set(h,'LineWidth',2);
set(h,'Marker','.','MarkerSize',20); % note you can specify multiple properties in one set() function call

A list of all the properties of a handle can be obtained by typing get(h).

MATLAB has a number of built-in handles, such as gca (“get current axes”):

set(gca,'LineWidth',2);
set(gca,'XTick',5904.62:0.01:5904.68,'YTickLabel',{});
set(gca,'FontSize',24);

☛ Look up what properties the gca handle has to find out how to change the color of the axes and labels, and change it to red.

Note that if you are making figures with multiple axes (perhaps by using subplot() or plotyy()) each axes gets its own handle. One axes is active at any given time, which it can be by clicking on it, or by explicitly doing something like axis(h).

There is also a handle for the whole figure, gcf:

set(gcf,'Color',[0 0 0]);

Finally, MATLAB has a special handle, '0', for some defaults. For instance, a very useful default to change is the font size, set(0,'DefaultAxesFontSize',18) because the default is usually too small to be readable once exported (see the next section on this). A good place for such a default change is in a shortcut that also contains your path, or perhaps even in MATLAB's startup.m.

### Exporting figures

Now that our figure looks nice, let's save it as an image file:

print(gcf,'-dpng','-r300','R042-2013-08-18-spikes.png');

Notice that the first argument of print is the current figure's handle. The other arguments specify that we want a PNG file, with 300dpi resolution, and the filename to write to. PNG format is a good choice for saving figures because it uses lossless compression, in contrast to JPEG images which use lossy compression and can have ugly artifacts as a result. Other useful save formats include -dill which saves .ai Illustrator files, and -deps which saves encapsulated PostScript; both of these are vector graphics formats.

☛ Look at the image file. Assuming you changed some of the figure and axis colors, you should find the colors in your image don't match those in MATLAB. Solve this problem by turning off the InvertHardCopy property of the figure, and save again. (MATLAB does this by default to facilitate printing images on white paper.)

You should get something like:

Obviously, I do not recommend formatting your rasterplots with this particular color scheme for publication!

### Putting it all together

Let's add some fancy plot features together:

%% restrict the data to interval of interest
this_iv = iv([5900 6000]);
S_r = restrict(S,this_iv);
csc_r = restrict(csc,this_iv);

%% plot spikes
SET_spkY = 0.4; % parameter to set spike height

figure; hold on;
for iC = 1:length(S_r.t)
if ~isempty(S_r.t{iC})
plot([S_r.t{iC} S_r.t{iC}],[iC-SET_spkY iC+SET_spkY],'k');
end
end % of cells
ylabel('neuron #');

SET_cscY = [-5 0];
plot(csc_r.tvec,rescale(csc_r.data,SET_cscY(1),SET_cscY(2)),'r');
set(gca,'YTickLabel',{});

%% add multi-unit activity in separate axes
ax1 = gca; ax1_pos = get(ax1,'Position');
ax2 = axes('Position',ax1_pos); % new axes with same position as first

cfg = []; cfg.tvec = csc.tvec; cfg.sigma = 0.1;
mua = getMUA(cfg,S); % obtain multi-unit activity

xr = get(ax1,'XLim');
mua_r = restrict(mua,xr(1),xr(2)); % only keep the data we need

axes(ax2); % set current axes to the second set, so what follows happens there
mua_hdl = plot(mua_r.tvec,mua_r.data,'Color',[0.7 0.7 0.7]);

set(gca,'YAxisLocation','right','Box','off','XTick',[],'Color','none','YColor',[0.5 0.5 0.5])
ylabel('multi-unit activity (spk/s)');
linkaxes([ax1 ax2],'x'); 

A few things to note here:

• Using the Position property of the current axes (ax1, containing the rasterplot) we created a new set of axes (ax2).
• In these new axes, we plotted the multi-unit activity, obtaining an output argument. Like gca and gcf this is a handle to graphics object; in this case, a handle to the multi-unit activity signal we just drew.
• We set the properties of the new multi-unit axes to have the y-axis on the right (plus a few other properties).
• linkaxes() was used to link the x-axis of both axes so that both update when zooming in.

☛ Use the 'XLim' axes property to zoom in to 5965 to 5969 s.

You should now be looking at the synchronous activation of a substantial number of neurons, reflected in the MUA peak, and associated with a high-frequency oscillation (~150-250Hz) in the LFP. These are the neurophysiological signatures of a “sharp wave-ripple complex” (SWR for short), events which are thought to contribute to the consolidation and retrieval of episodic memories.

It would be nice to be able to update what time window of the data we are looking at, without having to type these XLim commands. To do this we need to use a special Figure property, introduced in the next section.

The topics that follow in this section are optional, in the sense that later modules do not assume you know how to do these things. Feel free to use this section as you see fit, but do make sure you go through the next section (“Using existing visualization tools”)!

### Interactive figures and callback functions

Figure windows in MATLAB have many properties (here is the complete list). A particularly useful one is the KeyPressFcn property, which specifies a function to be called when a key is pressed while the figure is active.

For example, we can write a simple function that enables us to scroll left and right using the arrow keys, as follows:

function figscroll(src,event)

ax = get(src, 'CurrentAxes');
x_orig = get(ax, 'XLim');
x_step = (x_orig(1)-x_orig(2))/2;

switch event.Key
case 'leftarrow'
x_new = x_orig + x_step;
case 'rightarrow'
x_new = x_orig - x_step;
end

set(ax,'XLim',x_