~~DISCUSSION~~

===== Module 3: Visualizing neural data in MATLAB (short version) =====

Goals:

  * Load local field potential (LFP) and spike data from an example data set
  * Use the ''MultiRaster()'' function to plot and navigate through loaded raw data
  * Understand which file format to export figures to

Resources:

  * (reference) MATLAB documentation for ''plot'', [[ http://www.mathworks.com/help/matlab/ref/figure_props.html | Figure properties ]], [[ http://www.mathworks.com/help/matlab/ref/axes_props.html | Axis properties ]] and [[ http://www.mathworks.com/help/matlab/ref/line_props.html | Line Properties ]]
  * (reference) [[http://www.mathworks.com/help/matlab/creating_plots/figures-plots-and-graphs.html | MATLAB plotting gallery ]]
  * (reference) [[http://www.vandermeerlab.org/BIOL377_assignment3_2012.pdf | Place cell assignment]] from %%MvdM%%'s BIOL377 course (if you want a gentle intro to what is going on in the data you are looking at; can be done in 30 min)
  * (optional) Tufte, [[ http://www.edwardtufte.com/tufte/books_vdqi | The Visual Display of Quantitative Information ]]


==== Loading some data ====

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!). 

This data set was collected (by [[https://github.com/aacarey|Alyssa Carey]] for her Master's thesis work) 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.)


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

<code matlab>
%% cd to data folder -- replace this with yours
fc = 'D:\data\DataAnalysisTutorial\R042-2013-08-18';
cd(fc);

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

If you are using your own data that has a missing or bad .clusters file, do this instead:

<code matlab>
cfg = [];
cfg.useClustersFile = 0;
S = LoadSpikes(cfg)
</code>


This should give a **ts** data structure containing spike train data:

<code matlab>
LoadSpikes: Loading 67 files...

S = 

        t: {1x67 cell}
    label: {1x67 cell}
      cfg: [1x1 struct]
      usr: [1x1 struct]
</code>

The details of this are discussed in [[analysis:nsb2016: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:nsb2016: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:nsb2016:week2|Module 2]] a LFP **tsd** is defined by matching sets of sample timestamps (''csc.tvec'') and sampled data (''csc.data'').

==== Plotting basic spike rasters ====

The ''S'' variable that contains our spike times is a MATLAB struct that conforms to the ''ts'' ("timestamp") format introduced in [[analysis:nsb2016: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>

☛ 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.

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; we will do this below. Alternatively, there exists code that is much faster at plotting rasters (albeit more obscure; implemented in ''PlotSpikeRaster2()'' if you want to see an example).

==== Temporally restricting ts and tsd objects ====

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 a ''S_r'' variable with data restricted to between 5950 and 6050 s.

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

☛ 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 convenient data range (e.g. ''rescale(csc.data,-5,0)'').

You should see that the synchronous activation around time 5967 is indeed associated with a striking oscillation in the LFP:

{{ :analysis:nsb2014:m5_swr_example.png?600 |}}

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.

==== Exporting figures ====

A good way to save figures you may want to keep is:

<code matlab>
print(gcf,'-dpng','-r300','R042-2013-08-18-LFPsnippet.png');
</code>

Notice that the first argument of ''print'' is the current figure (''gcf''). 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). Other useful save formats are ''-dill'' which saves ''.ai'' Illustrator files, and ''-depsc'' which saves an embedded %%PostScript%% file (both vector formats).

☛ Critically evaluate the formatting of figure shown above. What are some missing elements you would want to see if this figure were to be used in a talk or paper? What changes could you make to improve legibility?

==== A fully featured plotting function ====

''MultiRaster()'' is a versatile and fast plotting function, with many different options and features.

The most basic usage is simply:

<code matlab>
MultiRaster([],S)
</code>

Note that the first argument in this case is the empty matrix ''[]''. For more features, this argument is used as a ''cfg'' variable with several options, for instance:

<code matlab>
plot_cfg = [];
plot_cfg.lfp = csc;
plot_cfg.spkColor = 'hsv';
MultiRaster(plot_cfg,S)
</code>

This particular data set also has some experimenter annotations associated with it, for example in the %%ExpKeys%% file:

<code matlab>
>> LoadExpKeys
>> ExpKeys

ExpKeys = 

                Behavior: 'MotivationalT'
         RestrictionType: 'water'
                 Session: 'standard'
                  Layout: 'foodLeft'
                Pedestal: 'R'
                nPellets: 5
             waterVolume: []
                 nTrials: 18
            forcedTrials: []
    nonConsumptionTrials: []
               badTrials: []
              pathlength: 318
                patharms: 369
           realTrackDims: [139 185]
                convFact: [2.9176 2.3794]
             TimeOnTrack: 3240
            TimeOffTrack: 5645
               prerecord: [2.1266e+03 3.2141e+03]
                    task: [3.2387e+03 5.6452e+03]
              postrecord: [5.6564e+03 6.5635e+03]
                 goodSWR: {1x2 cell}
               goodTheta: {'R042-2013-08-18-CSC07a.ncs'}

</code>

The meaning of these fields is explained in more detail in the full [[task description]], but for now we will just use the start and end times for the 'prerecord' segment of the data to highlight it visually:

<code matlab>
plot_cfg = [];
plot_cfg.evt = iv(ExpKeys.prerecord);
plot_cfg.evtColor = [0 0 0.5];
MultiRaster(plot_cfg,S)
</code>

''MultiRaster()'' is interactive, allowing the user to scroll, zoom, move from one event to another, et cetera. To see the full list of keybindings, type ''help navigate''.