Differences

This shows you the differences between two versions of the page.

--- analysis:nsb2014:week5 [2014/07/20 12:48]
mvdm [Module 5: Visualizing neural data in MATLAB]
+++ analysis:nsb2014:week5 [2018/07/07 10:19] (current)
@@ Line 34: / Line 34: @@
 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:
+==== Loading some data ====
-==== Data pre-processing overview ====
+First, 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!).
-First, make sure you have the folder ''R042-2013-08-18'' from the lab database, 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. 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.)
@@ Line 49: / Line 47: @@
 <code matlab>
+%% cd to data folder -- replace this with yours
+fd = 'D:\data\DataAnalysisTutorial\R042-2013-08-18';
+cd(fc);
 %% load the data (note, may need to unzip position data first)
-fc = FindFiles('*.t');
+cfg = [];
-S = LoadSpikes(fc);
+S = LoadSpikes(cfg)
+</code>
-[csc,csc_info] = LoadCSC('R042-2013-08-18-CSC03a.ncs');
+If you are using your own data that has a missing or bad .clusters file, do this:
-[Timestamps, X, Y, Angles, Targets, Points, Header] = Nlx2MatVT('VT1.nvt', [1 1 1 1 1 1], 1, 1, [] );
+<code matlab>
+cfg = [];
+cfg.useClustersFile = 0;
+S = LoadSpikes(cfg)
 </code>
-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.
+This should give a data structure containing spike train data:
-  * ''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 ====
+<code matlab>
+LoadSpikes: Loading 67 files...
-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.
+S =
-''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.
+        t: {1x67 cell}
+    label: {1x67 cell}
+      cfg: [1x1 struct]
+      usr: [1x1 struct]
+</code>
-''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()''.
+The details of this are discussed in [[analysis:nsb2014:week2|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 [[analysis:nsb2014:week4|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.
+Let's also load a LFP:
+<code matlab>
+cfg = [];
+cfg.fc = {'R042-2013-08-18-CSC03a.ncs'};
+csc = LoadCSC(cfg);
+</code>
+This gives:
+<code matlab>
+LoadCSC: Loading 1 files...
+LoadCSC: R042-2013-08-18-CSC03a.ncs 0/17193 bad blocks found (0.00%).
+>> csc
+csc =
+     tvec: [8802816x1 double]
+     data: [1x8802816 double]
+    label: {'R042-2013-08-18-CSC03a.ncs'}
+      cfg: [1x1 struct]
+</code>
+''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 in [[analysis:nsb2014:week2|Module 3]] a LFP is defined by matching sets of sample timestamps (''csc.tvec'') and sampled data (''csc.data'').
+==== ts objects: basic methods ====
+The ''S'' variable that contains our spike times is a MATLAB struct that conforms to the ''ts'' ("timestamp") format introduced in [[analysis:nsb2014:week2|Module 2]]. Its ''.t'' field contains spike trains as lists of numbers, corresponding to the timestamps of the spikes. So, to plot a raster with the spikes of the first cell in our data set:
+<code matlab>
+nCell = 1;
+SET_SpikeHeight = 0.4;
+plot([S.t{nCell} S.t{nCell}],[nCell-SET_SpikeHeight nCell+SET_SpikeHeight],'k')
+set(gca,'YLim',[-10 10]);
+</code>
-☛ Verify that your ''csc'' object indeed has timestamps and data of the same length.
+☛ Extend the example code above to make a rasterplot containing the spikes for all cells in ''S''. To facilitate re-use later, wrap your code in a function, %%PlotSpikeRaster%%, that plots the raster for any ''ts'' input.
-By default the timestamps are given in seconds (s).
+You should find that it takes a significant amount of time to plot that many spikes. In practice we often want to plot spikes for specific time intervals.
-A useful method that works both on ''ts'' and ''tsd'' objects is ''Restrict()''.
+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 cell, use it to create a ''csc_short'' variable with data restricted to between 5950 and 6050 s.
+☛ Use MATLAB's ''help'' function to look up the usage information for ''restrict()''. In a new editor cell, use this function to create a ''S_r'' 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.
+Calling your %%PlotSpikeRaster%% function on this restricted data should be much faster. Now that we are a bit more zoomed in, some interesting features of the rasterplot become visible. For instance, at time 5967, many cells are synchronously active. Could this be a sharp wave-ripple event (SWR)?
-☛ Look at the ''Timestamps'' variable and notice that it's unlikely to be in seconds.
+☛ Restrict the csc to the same time (5950 to 6050) and plot it in the same figure. You can use the ''rescale()'' function to fit the csc into a convienient data range (e.g. ''rescale(csc.data,-5,0)'').
-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()''.
+You should see that the synchronous activation around time 5967 is indeed associated with a striking oscillation in the LFP:
-If you get an error, try creating the tsd with a transposed data argument (''X''')
+{{ :analysis:nsb2014:m5_swr_example.png?600 |}}
-==== Plotting with handles ====
+This high-frequency osillation (the "ripple") superimposed on a slower oscillation (the "sharp wave") is the signature associated with "replay" events in the hippocampus. We will explore this phenomenon in detail in later modules.
+==== Figure properties, handles, and exporting ====
-☛ 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.
+(Remember to execute the commands below from your ''sandbox.m'' file, not from the command line! That way you can more easily retrieve what you did later.)
-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
+A 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
 <code matlab>
@@ Line 100: / Line 145: @@
 ''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).
+☛ 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.
@@ Line 116: / Line 161: @@
 ☛ 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.
+Let's add some moderately fancy features to our plot. Starting from a currently selected rasterplot, do the following:
 <code matlab>
-hold on; box off;
+%% add multi-unit activity to rasterplot
+ax1 = gca; ax1_pos = get(ax1,'Position');
+ax2 = axes('Position',ax1_pos); % new axes with same position as first
-csc_mean = nanmean(Data(csc));
+cfg = []; cfg.tvec = csc.tvec; cfg.sigma = 0.1;
-xr = get(gca,'XLim');
+mua = getMUA(cfg,S); % obtain multi-unit activity
-mean_hdl = plot(xr,[csc_mean csc_mean]);
+xr = get(ax1,'XLim');
+mua_r = restrict(mua,xr(1),xr(2)); % only keep the data we need
+axes(ax2);
+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')
+linkaxes([ax1 ax2],'x');
 </code>
 A few things to note here:
-  * ''hold on'' prevents our newly plotted object from erasing whatever was there before.
+  * Using the ''Position'' property of the current axes (ax1, containing the rasterplot) we created a new set of axes (ax2).
-  * ''box off'' speaks for itself.
+  * 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 used the ''XRange'' property of the current axes to get the smallest and largest x values in the plot.
+  * We set the properties of the new multi-unit axes to have the y-axis on the right (plus a few other properties).
-  * ''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.
+  * ''linkaxes()'' was used to link the x-axis of both axes so that both update when zooming in. Try it!
-☛ Change the properties of ''mean_hdl'' so that the line is red, dashed, and of width 2.
+☛ Change the properties of ''mua_hdl'' so that the line is red, dashed, and of width 2.
 seconds of LFP is a lot to display.
@@ Line 139: / Line 193: @@
 ☛ 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.
+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.
 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.
@@ Line 195: / Line 249: @@
 ☛ 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 ====
+A different way of flexibly specifying function inputs is through the use of a ''cfg'' (configuration) input. An example of how this works can be seen in the ''getMUA()'' function used above.
-☛ Implement a function ''neuroplot()'' as follows:
+==== Exercise ====
+☛ Extend your %%PlotSpikeRaster%% function as follows, using either the varargin or cfg strategy, and your knowledge of object handles.
 <code matlab>
-function neuroplot(spikes,csc,varargin)
+% required inputs:
-% function neuroplot(spikes,csc,varargin)
 %
-% inputs:
+% spikes: ts object with spike train
 %
-% spikes: {nCells x 1} cell array of ts objects (spike train)
+% optional inputs:
-% csc: {nCSCs x 1} cell array of tsd objects (LFPs)
 %
-% varargins:
+% csc: tsd object with LFP data (bonus points for accepting more than one gracefully)
 %
-% cscColor: [nCSC x 3] RGB values to plot CSCs, default []
+% cscColor: [nCSC x 3] RGB values to plot CSCs, defaults to red
-% spikeColor: [nCells x 3] RGB values to plot spikes, default []
+% spikeColor: [nCells x 3] RGB values to plot spikes, defaults to black
 %
-% evt: [nEvents x 1] event times to plot, default []
+% evt: ts object with event times to plot as vertical lines, defaults to none
-% evtColor: [nEvents x 3] RGB values to plot events, default []
 %
-% interactiveMode: boolean to enable/disable arrow key navigation
+% interactiveMode: boolean to enable/disable arrow key navigation of time axis
 </code>
-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.
+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. 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 (such as times when the animal received reward, or brain stimulation was triggered); 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.
-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.
+☛ When you are satisfied with your work, save a copy in your project folder and push it to GitHub. Verify that your version is visible from the web.
 === Hints ===

User Tools

Site Tools

Differences

Page Tools