User Tools

Site Tools


analysis:course:week2

Visualizing neural data in MATLAB

Goals:

  • Load LFP and spike data from an example data set
  • Master basic and intermediate MATLAB plotting: plot, axis and figure properties
  • Use interactive figure properties for dynamic data viewing
  • Learn to use functions with variable numbers of input arguments (varargin) effecively
  • Understand which file format to export figures to

Deliverables:

  • A nice rasterplot function with optional arguments to add events, color, and interactive navigation.

Resources:

Introductory remarks

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 touching up and combining them in Illustrator. We will do the MATLAB part here and aim to produce something similar to panel Aa.

An important principle in graphical design for scientific communication is that you want to only include those elements which contribute to the messages(s) you want the graphic to communicate. Remove and simplify as much as possible: this will make your figures more readable. The Tufte book 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 your figure legends big enough to be readable! :-)

Step-by-step

First, some housekeeping.

Use the shortcut button created in the previous module to restore MATLAB's default path, add the vandermeerlab codebase to it, and move to your BIOL680 folder.

Create a folder with today's date and cd to it. This will be your “scratch” working directory with stuff that you don't necessarily intend to push to GitHub (although you are welcome to if it's easier or you would like me to be able to see what you did). At the end of this module you will have something to push, but this will go in the shared folder under BIOL680 once you are satisfied with your work.

Grab the folder R042-2013-08-18 from the lab database, making sure you login with BIOL680 credentials so we all have the same version. As before, place the data in a sensible local location. 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. 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.)

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

%% load the data (note, may need to unzip position data first)
fc = FindFiles('*.t');
S = LoadSpikes(fc);
 
[csc,csc_info] = LoadCSC('R042-2013-08-18-CSC03a.ncs');
 
[Timestamps, X, Y, Angles, Targets, Points, Header] = Nlx2MatVT('VT1.nvt', [1 1 1 1 1 1], 1, 1, [] );

You now have a number of variables in your workspace. The most important ones are:

  • 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. There is also csc_info, a struct with a some properties of the csc you just loaded. We will explore this in detail in the next module. Notice that csc is a tsd object; more on this in the next section.
  • S is a cell array of ts objects (discussed in the next section). Each cell of S contains a ts object, a spike train 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.
  • X and Y contain the location of the rat as tracked by an overhead camera. Each element has a corresponding time in Timestamps.

ts and tsd objects: basic methods

MATLAB supports basic “Object-oriented programming” functionality. If you are not familiar with this, it is sufficient for now to know that objects, such as those of the ts and tsd class you have in your workspace, have specific functions associated with them. These class-specific functions or methods generally do not work on objects of other classes, or variables that are not objects. ts and tsd are classes not included with MATLAB; they are part of the lab codebase.

ts stands for “timestamp”. The data in a ts object, which you can access with Data() is simply a list of numbers. For our spike trains here these numbers indicate the times at which that neuron fired a spike.

tsd stands for “timestamped data”. For a ts object, the timestamps fully describe the data. Often however we want to describe data where timestamps are associated with a quantity such as a voltage or a position: timestamped data. For our csc, we can access the data with Data() and the matching timestamps with Range().

☛ Verify that your csc object indeed has timestamps and data of the same length.

By default the timestamps are given in seconds (s).

A useful method that works both on ts and tsd objects is Restrict().

☛ Use MATLAB's help function to look up the usage information for Restrict. In a new cell, use it to create a csc_short variable with data restricted to between 5950 and 6050 s.

Restrict() ensures that timestamps and data never get out of sync, a useful feature during analysis. Let's create a tsd object for our X and Y position data.

☛ Look at the Timestamps variable and notice that it's unlikely to be in seconds.

In fact, it is in Neuralynx's internal clock units. After converting to seconds (multiply by 10^-6), create a tsd for X and Y by using the tsd() function for each. Verify that it worked by creating a X_short variable over the same range as csc_short by using Restrict().

If you get an error, try creating the tsd with a transposed data argument (X')

Plotting with handles

☛ Plot the LFP (time against voltage) for the segment between 5950 and 6050 s. Remember to do this using a cell in your sandbox.m file so that it's easy to make changes later and you keep a record of what you did.

Any MATLAB figure is actually a collection of objects, all with their own properties which you can access and edit. For instance, to change the figure background to black, do

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

gcf is MATLAB's name for “get current figure”. The above command set its 'Color' property to [0 0 0] which is black (zero red, zero green, zero blue in the RGB world).

☛ Set the background of the current /axes/ to black (gca is MATLAB's name for the current axes object; notice that the axes background is currently white, which looks ugly).

Notice that now you can't see the axis lines and labels anymore, because they too became black when you changed the color of the axes. You want to set their color to be white, but you might not know the name of the relevant property.

☛ Use get(gca) to get a list of all the axes properties. Find the properties that correspond to the axes color and set them to white.

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

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

Notice that the first argument of print is the current figure. The other arguments specify that we want a PNG file, with 300dpi resolution. PNG format is a good choice for saving figures because it uses lossless compression, JPG images, which use lossy compression, can have ugly artifacts (but more colors). Another useful save format is -dill which saves .ai Illustrator files.

☛ Look at the image file. The colors are all wrong! 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.)

Let's add some more features to our plot.

hold on; box off;
 
csc_mean = nanmean(Data(csc));
xr = get(gca,'XLim');
 
mean_hdl = plot(xr,[csc_mean csc_mean]);

A few things to note here:

  • hold on prevents our newly plotted object from erasing whatever was there before.
  • box off speaks for itself.
  • We used the XRange property of the current axes to get the smallest and largest x values in the plot.
  • plot() was called with an output argument. Like gca and gcf this is a handle to graphics object; in this case, a handle to the line we just drew.

☛ Change the properties of mean_hdl so that the line is red, dashed, and of width 2.

100 seconds of LFP is a lot to display.

☛ Use the 'XLim' axes property to zoom in to 5989 to 5990 s; while you are at it, increase the font size to 24. (You can change multiple properties with the same command, just keep adding property-value pairs.)

It would be nice to be able to update the figure without having to type these XLim commands. To do this we need to use a special Figure property, introduced in the next section.

One final point of interest: different graphics objects are often (hierarchically) related. For instance, our figure currently contains one set of axes. We can obtain the axes handle directly, axes_handle = gca;, but also indirectly, axes_handle = get(gcf,'Children');. This works because currently, the axes are the only child of the figure.

Figure callback functions

First, another basic MATLAB skill: anonymous functions. Consider:

sqr_fn = @(x) x.^2;
sqr_fn(2)

sqr_fn defines a function handle to a so-called anonymous function, i.e. a function that is not in .m file form. The above example basically says “sqr_fn is a function of x, and should return x squared.” (if you don't understand the . notation, it is important to look this up!)

Look at the KeyPressFcn figure property. This property specifies the function to call whenever a key is pressed with the figure active..

☛ Implement the example shown in the page above. Notice the properties of the evt object passed to the keypress function, and that this function also gets passed the handle of the figure that generated the keypress event.

The keypress function thus has access to everything we need to implement some useful navigation functionality.

☛ Create a keypress function that moves the x-limits of the current figure's axes to the left or right by 50% (i.e. if the current axes are from 1 to 2s, then pressing the right arrow key takes it to 1.5-2.5 s. It will be easier to do this with a .m file given the length of what you need to do.

Hint: model your .m file on the template in the example given in the KeyPressFcn section of the Figure properties page, and define it thus: f_hdl = figure('KeyPressFcn',@your_function); if your function is called your_function.m. Please keep in mind of the difference between ”==” and “strcmp” operators. The former treats two strings as vectors and does letter-wise comparison. Thus the two strings ought to be of the same dimension. The latter operator does lexical comparison and returns “1” if the two strings are identical, “0” otherwise.

Varargins

A standard function definition specifies the exact number of input and output arguments the function expects: function y = sqrt(x) has one input and one output argument. However, it is often useful to have a variable number of input arguments, so that you can override defaults or specify additional options as needed without complicating simple function calls. For instance, you want to be able to make quick plots with plot(x) but specify extra stuff when needed (plot(x,'LineWidth',2)). This is accomplished with the special argument varargin.

Here is the idea:

function test_fun(a,varargin)
 
b = 2; % set defaults for b and c
c = 3;
 
extract_varargin; % override b and c if specified
 
fprintf('a is %d, b is %d, c is %d\n',a,b,c);

Test it:

>> test_fun(1)
a is 1, b is 2, c is 3
>> test_fun(1,'b',1)
a is 1, b is 1, c is 3
>> 

Note how specifying a “key-value” pair ('b',1) overwrote the default for b. extract_varargin basically goes through the varargins and assigns each value to the corresponding key.

☛ Verify your understanding by calling test_fun such that b and c are both set to 0. Also note what happens when you do test_fun(1,'B',1)!

Assignment

☛ Implement a function neuroplot() as follows:

function neuroplot(spikes,csc,varargin)
% function neuroplot(spikes,csc,varargin)
%
% inputs:
%
% spikes: {nCells x 1} cell array of ts objects (spike train)
% csc: {nCSCs x 1} cell array of tsd objects (LFPs)
%
% varargins:
%
% cscColor: [nCSC x 3] RGB values to plot CSCs, default []
% spikeColor: [nCells x 3] RGB values to plot spikes, default []
%
% evt: [nEvents x 1] event times to plot, default []
% evtColor: [nEvents x 3] RGB values to plot events, default []
%
% interactiveMode: boolean to enable/disable arrow key navigation

This function should produce a rasterplot of spikes – one neuron per row – with one or multiple LFPs plotted above it, similar to panel Aa of the Dragoi and Tonegawa figure above. The values on the y-axis are not important for now, so you can feel free to hide this axis and rescale everything to look good. (However, if you want to include a scalebar for the LFPs, go ahead!) Set the default x range to 2 seconds.

You can pick a default color to plot everything, but there should be optional arguments to have the user specify a color for each neuron and LFP.

Also, we want to be able to overlay events of interest on this plot; these can be simple vertical lines that extend all the way through the plot. You can make up a few event times yourself for testing.

Finally, for data exploration it is helpful to be able to skip through the data with the arrow keys; set it up so that the left and right arrows can be used to scroll left and right. For extra awesomeness, add zoom functionality that scales the LFPs when pressing the up and down arrows.

☛ When you are satisfied with your work, save a copy in your shared folder and push it to GitHub. Verify that your version is visible from the web.

Hints

Everything you need to know to make this work has been introduced in this module: beyond this (and basic MATLAB) no special functions or tricks are needed. A few suggestions:

  • Start simple. Instead of trying to implement everything at once, take it one step at a time. First, plot a single line (corresponding to one spike). Then, plot a spike train for one neuron, and keep building.
  • Don't test on the full data set. Things will be faster if you use restricted data for testing.
  • Plot everything in a single set of axes, i.e. without using subplot() or axes(); this means you have to rescale and add a baseline offset (a.k.a. DC component) to the LFP signals so that they don't overlap.
  • Develop in the sandbox file so that you can use cell mode and have easy access to the workspace. Only when everything is basically working, move to putting everything inside a function.

Discussion

Syaheed Jabar, 2013/09/19 11:45

Just a note/query:

Is there an issue with trying to implement the KeyPressFcn function the way the Matlab documentation suggests (using a nested anonymous function), while also using the varargin function so the user can specify additional arguments?

Seems like using nested functions requires the ‘end’ statement to the function, which seemingly creates a static workspace, which then seemingly disables the ability to use the variable number of arguments that varargin supposedly implements. Not using nested functions (putting the previously nested function in another m-file seems to side-step the issue, however.

Just found it a bit curious, although it could just be that I am mistaken …

Matt van der Meer, 2013/09/20 11:59, 2013/09/20 12:05

I don't see how the use of varargin would interact with any KeyPressFcn unless you would like the KeyPressFcn to accept a variable number of input arguments. Note that MATLAB only allows anonymous function definitions, not other function definitions, in scripts

Derek Zhi, 2013/09/24 13:16

With regards to the Figure callback functions, it would be nice to distinguish the difference between ”==” and “strcmp” operators. The former treats two strings as vectors and does letter-wise comparison. Thus the two strings ought to be of the same dimension. The latter operator does lexical comparison and returns “1” if the two strings are identical, “0” otherwise.

Matt van der Meer, 2013/09/24 13:20

Good point – the example given in the MATLAB documentation is actually a bit sloppy in using ==! Go ahead and note this in the main document text above if you wouldn't mind.

Eric Carmichael, 2013/09/24 19:38

Has anyone identified any specific evt times in the data set? Am I right in thinking that these times should lead to the non-continuous sections seen in the Dragoi figure? So far I only have a continuous raster and LFP with a specific x limit to be displayed but that does not lead to these sections.

Matt van der Meer, 2013/09/24 21:02

You don't need to produce a multipanel figure as in the Dragoi example – just one panel (so one “box” in the Aa part of the figure) is sufficient. You aren't given any event times, so you can just make some up. For instance, something like every 10 seconds is fine. It is just to demonstrate that the function could accept event inputs and plot them in a useful way. Sorry for any confusion on this point!

Matt, 2021/09/25 07:51

A leisurely walking tour through the historical center of the capital is a special kind of leisure, equally interesting for guests and residents of the city. The classic route includes sightseeing of iconic sights and architectural masterpieces, contemplation of panoramic views and features of Moscow. You will see how people work and relax, take memorable photos, you will be able to enjoy the color of the city on its streets and parks.

Enter your comment. Wiki syntax is allowed:
 
analysis/course/week2.txt · Last modified: 2018/07/07 10:19 (external edit)