Differences

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

--- analysis:nsb2014:week8 [2014/07/23 11:59]
mvdm
+++ analysis:nsb2014:week8 [2022/06/29 19:01] (current)
mvdm
@@ Line 1: / Line 1: @@
-~~DISCUSSION~~
-Under construction!
 ===== Time-frequency analysis: spectrograms =====
@@ Line 49: / Line 45: @@
 Let's give it a spin:
-FIXME
 <code matlab>
 %% load the data
 % remember to cd to the correct folder here, may need to get this file from the lab database
-fname = 'R016-2012-10-03-CSC04a.Ncs';
+cfg = []; cfg.fc = {'R016-2012-10-03-CSC04a.Ncs'};
-csc = myLoadCSC(fname);
+csc = LoadCSC(cfg);
-hdr = getHeader(csc);
-Fs = hdr.SamplingFrequency;
+Fs = csc.cfg.hdr{1}.SamplingFrequency;
 %% restrict data
-cscR = Restrict(csc,3282,3286); % if you don't have this, implement it (it's one line of code!)
+cscR = restrict(csc,3282,3286); % if you don't have this, implement it (it's one line of code!)
-plot(cscR); % note there are various gamma oscillations present, as well as a large negative-going transient
+plot(cscR.tvec,cscR.data); % note there are various gamma oscillations present, as well as a large negative-going transient
 %% construct and plot the spectrogram
-[S,F,T,P] = spectrogram(Data(cscR),hanning(512),256,1:200,Fs);
+[S,F,T,P] = spectrogram(cscR.data,hanning(512),256,1:200,Fs);
 imagesc(T,F,10*log10(P)); % converting to dB as usual
 set(gca,'FontSize',20);
@@ Line 94: / Line 87: @@
 <code matlab>
->> tvec = Range(csc_postRec);
+>> cscR.tvec(1)
->> tvec(1)
 ans =
@@ Line 102: / Line 94: @@
 </code>
-In fact, the time of the first bin returned by ''spectrogram()'' is the //centre of the first window for which there is enough data//. The data passed to ''spectrogram()'' is assumed to start at time 0. This makes sense when we remember that the input arguments do not contain timestamps (only voltage samples) so this function has no way of knowing the "experiment time" (''Range(csc)'').
+In fact, the time of the first bin returned by ''spectrogram()'' is the //center of the first window for which there is enough data//. The data passed to ''spectrogram()'' is assumed to start at time 0. This makes sense when we remember that the input arguments do not contain timestamps (only voltage samples) so this function has no way of knowing the "experiment time" (''csc.tvec'').
-Thus, the time of the first spectrogram bin (''0.128'') is because this is the centre of the first 512-sample window for which there is enough data (''256/2000'' is ''0.128''), and the first sample of the data is assigned time 0.
+Thus, the time of the first spectrogram bin (''0.128'') is because this is the center of the first 512-sample window for which there is enough data (''256/2000'' is ''0.128''), and the first sample of the data is assigned time 0.
 === Aligning the spectrogram with a LFP ===
@@ Line 111: / Line 103: @@
 <code matlab>
+%%
 hold on;
-tvec = Range(cscR);
-data = Data(cscR);
 lfp_minmax = 25; lfp_cent = 125; % range and mean of LFP plotting
-tvec0 = tvec - tvec(1); % align LFP with spectrogram
+tvec0 = cscR.tvec - cscR.tvec(1); % align LFP with spectrogram
-data = rescale(data,-lfp_minmax,lfp_minmax); data = data+lfp_cent;
+data = rescale(cscR.data,-lfp_minmax,lfp_minmax); data = data+lfp_cent;
 lfp_h = plot(tvec0,data,'k');
 </code>
@@ Line 177: / Line 168: @@
 <code matlab>
-cscR = Restrict(csc,3300,3340); % section of data with a gap
+cscR = restrict(csc,3300,3340); % if you don't have this, implement it (it's one line of code!)
-[S,F,T,P] = spectrogram(Data(cscR),rectwin(256),128,1:200,Fs);
+[S,F,T,P] = spectrogram(cscR.data,rectwin(256),128,1:200,Fs);
+imagesc(T,F,10*log10(P)); % converting to dB as usual
-imagesc(T,F,10*log10(P));
 set(gca,'FontSize',20);
 axis xy; xlabel('time (s)'); ylabel('Frequency (Hz)');
 hold on;
-tvec = Range(cscR); data = Data(cscR);
 lfp_minmax = 25; lfp_cent = 125; % range and mean of LFP plotting
-tvec0 = tvec - tvec(1); % align LFP with spectrogram
+tvec0 = cscR.tvec - cscR.tvec(1); % align LFP with spectrogram
-data = rescale(data,-lfp_minmax,lfp_minmax); data = data+lfp_cent;
+data = rescale(cscR.data,-lfp_minmax,lfp_minmax); data = data+lfp_cent;
 lfp_h = plot(tvec0,data,'k');
 xlim([tvec0(1) tvec0(end)]);
@@ Line 216: / Line 205: @@
 === Processing Neuralynx events files (.nev) ===
-Apart from neural data stored in ''.ncs'' and ''.ntt'' files, the Neuralynx system also stores //events// in a file called ''Events.Nev''. Let's load it to see what is inside:
+Apart from neural data stored in ''.ncs'' and ''.ntt'' files, the Neuralynx system also stores //events// in a file called ''Events.Nev'' (see [[analysis:nsb2014:week2|Module 2]] for details). Here is an example usage for the specific task in this data:
 <code matlab>
-fn = FindFile('*Events.nev');
+cfg = [];
-[EVTimeStamps, EventIDs, TTLs, EVExtras, EventStrings, EVHeader] = Nlx2MatEV(fn,[1 1 1 1 1],1,1,[]);
+cfg.eventList = {'Feeder 0 nosepoke','Feeder 1 nosepoke', ...
+    '1 pellet cue','3 pellet cue','5 pellet cue', ...
+    '1 pellet dispensed','3 pellet dispensed','5 pellet dispensed'};
+cfg.eventLabel = {'n0','n1', ...
+    'c1','c3','c5', ...
+    'd1','d3','d5'};
+evt = LoadEvents(cfg);
 </code>
-If you look at the ''EventStrings'' variable, you will see something like:
+This result is the following:
 <code matlab>
-    'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0004).'
+>> evt
-    '1 or 5 pellet cue'
-    'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0006).'
-    'Feeder 1 nosepoke'
-    'TTL Output on AcqSystem1_0 board 0 port 0 value (0x0002).'
-    '1 pellet dispensed'
-    'TTL Output on AcqSystem1_0 board 0 port 0 value (0x0000).'
-    'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0004).'
-    'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0006).'
-    'TTL Input on AcqSystem1_0 board 0 port 1 value (0x0004).'
-</code>
-Each //cell// of ''EventStrings'' contains a string that corresponds to an event. Events are stored based on a number of different triggers, including:
+evt =
-  * Any change that is detected on the input/output (I/O) ports of the system; these ports interface with components of the experimental setup such as reward pellet dispensers, light cues, levers, photosensors, and others
+        t: {1x8 cell}
-  * System events, such as "Start Recording", "Data lost", et cetera
+    label: {'n0'  'n1'  'c1'  'c3'  'c5'  'd1'  'd3'  'd5'}
-  * Manual keyboard input into the Cheetah software
+      cfg: [1x1 struct]
-  * Events sent from MATLAB
-In the above example, the "TTL Input..." and "TTL Output..." entries are all automatically generated because an input was detected (e.g. a photobeam break or lever press) or an output was sent (e.g. dispense a reward pellet). Messages such as ''Feeder 1 nosepoke'' were sent from the MATLAB script that controls the experiment. If you look through the ''EventStrings'' you will also find some system messages.
+>> evt.t
-The ''EventStrings'' variable has the same size as ''EVTimeStamps'' which contains the corresponding times for each event. This way, we can line up the events with the neural data.
+ans =
-=== An example getEvents.m file ===
+  Columns 1 through 4
-Typically we are interested in a specific set of events, such as all times of reward receipt, and want to extract these from the Events file. For the current task, the ''getEvents_value.m'' function extracts a number of events of interest:
+    [1x58 double]    [1x57 double]    [1x35 double]    [1x19 double]
-<code matlab>
+  Columns 5 through 8
->> evt = getEvents_value
-evt =
+    [1x34 double]    [1x40 double]    [1x25 double]    [1x42 double]
-    c1: [1x35 double]
-    c3: [1x19 double]
-    c5: [1x34 double]
-    d1: [1x40 double]
-    d3: [1x25 double]
-    d5: [1x42 double]
-    n0: [1x58 double]
-    n1: [1x57 double]
 </code>
-As you can see, this function returns a struct with a number of fields. ''c_x'' are the times of three different audio cues, indicating the availability of 1, 3 and 5 food pellets; ''d_x'' are the times at which the corresponding number of pellets were dispensed, and ''n_x'' are the times at which the rat nosepoked in receptacle 0 (on the left end of a linear track) or 1 (on the right).
+Thus, for each item in ''cfg.eventList'' we extract the corresponding event timestamps, placing them in an output structure field with a name specified in ''cfg.eventLabel''.
+The events returned by ''LoadEvents()'' above are already present in the Events file; this function merely presents them in a more usable format. However, for many tasks, some more processing is needed to obtain event times of interest, for instance "only those nosepokes that follow presentation of a cue predicting 5 pellets, and for which the rat does not un-poke for at least one second". Creating appropriate "trial functions" that return usable event times is unglamorous, but a key analysis step for any data set. (Remember: "Garbage in, garbage out"!)
-The events returned by ''getEvents_value()'' are already present in the Events file; this function merely presents them in a more usable format. However, for many tasks, some more processing is needed to obtain event times of interest, for instance "only those nosepokes that follow presentation of a cue predicting 5 pellets, and for which the rat does not un-poke for at least one second". Creating appropriate "trial functions" that return usable event times is unglamorous, but a key analysis step for any data set. (Remember: "Garbage in, garbage out"!)
-In any case, armed with some event times, we are now in a position to ask questions of the sort, "does reward delivery elicit changes in the ventral striatal LFP?" and "is the ventral striatal LFP sensitive to reward-predictive cues?" To do so, we will use the FieldTrip toolbox, which massively expands MATLAB's suite of built-in functions for spectral analysis.
+In any case, armed with some event times, we are now in a position to ask questions of the sort, "does reward delivery elicit changes in the ventral striatal LFP?" and "is the ventral striatal LFP sensitive to reward-predictive cues?" To do so, we will use the %%FieldTrip%% toolbox, which massively expands MATLAB's suite of built-in functions for spectral analysis.
 ==== Event-triggered spectrograms using FieldTrip ====
@@ Line 282: / Line 259: @@
 === FieldTrip setup ===
-☛ Download FieldTrip from the lab database using the BIOL680 account. For this course, do not use the official version from the FieldTrip website: we currently run a custom version modified to allow the processing of Neuralynx data. Unzip the archive in a suitable code subfolder -- do not mix code with your data!
+☛ If you haven't already done so, clone the lab's [[https://github.com/mvdm/fieldtrip|FieldTrip repository on GitHub]]. We currently run a custom version modified to allow the processing of Neuralynx data, so the official %%FieldTrip%% release will not work! See [[analysis:nsb2014:week1|Module 1]] for more details on setting this up.
-☛ Next, create a new MATLAB shortcut that is a copy of your existing course shortcut, but with the addition of the fieldtrip folder (and all its subfolders) to the path in addition to the GitHub codebase. I have named this shortcut ''mvdmlab-ft''.
-Click this shortcut to add FieldTrip to your path. If you can type ''ft_defaults'' and not get errors, the installation is successful.
+☛ Also make sure you have a MATLAB shortcut to create a path that includes the nsb-2014 and fieldtrip code (but not %%MClust%%!). Run this shortcut before proceeding.
 === Loading Neuralynx data into FieldTrip ===
@@ Line 321: / Line 295: @@
 === Creating trial data in FieldTrip ===
-To convert this single "trial" containing the entire recording into proper trials surrounding events of interest, FieldTrip uses a "trialfun": a function that returns the times of the events of interest. These times, along with some parameters like how much time around the events to include, are then passed to the ''ft_redefinetrial()'' function, which cuts up the data into trials:
+To convert this single "trial" containing the entire recording into proper trials surrounding events of interest, %%FieldTrip%% requires a ''trl'' variable that specifies the start and end indices (into the data) which are then passed to the ''ft_redefinetrial()'' function, which cuts up the data into trials:
 <code matlab>
 cfg = [];
-cfg.trialfun = 'ft_trialfun_lineartracktone2';
+cfg.t = cat(2,getd(evt,'n0'),getd(evt,'n1'));
-cfg.trialdef.hdr = data.hdr;
+cfg.mode = 'nlx';
-cfg.trialdef.pre = 2.5;
+cfg.hdr = data.hdr;
-cfg.trialdef.post = 5;
+cfg.twin = [-1 4];
-cfg.trialdef.eventtype = 'nosepoke'; % could be 'nosepoke', 'reward', 'cue'
+trl = ft_maketrl(cfg);
-cfg.trialdef.location = 'both'; % could be 'left', 'right', 'both'
-cfg.trialdef.block = 'both'; % could be 'value', 'risk'
-cfg.trialdef.cue = {'c1','c3','c5'}; % cell array with choice of elements {'c1','c3','c5','lo','hi'}
-[trl, event] = ft_trialfun_lineartracktone2(cfg);
+cfg = [];
 cfg.trl = trl;
+data_trl = ft_redefinetrial(cfg,data);
-data_trl = ft_redefinetrial(cfg,data);
 </code>
-The above reorganizes the data into 7-second windows (from -2 to 5 seconds) around nosepoke times following the 1- 3- and 5-pellet cues. The function ''ft_trialfun_lineartracktone2()'' is a more powerful version of ''getEvents_value()'' we encountered above: it can accept a variety of requests for particular event times. For instance, if we only wanted nosepokes on the left side of the maze that followed the 3-pellet cue, we can specify this by changing 'both' to 'left' and only including 'c3'.
+The above reorganizes the data into 4-second windows (from -1 to 4 seconds) around nosepoke times (''evt.n0'' and ''evt.n1''). ''ft_maketrl()'' converts the event times from ''LoadEvents()'' into the data indices that ''ft_redefinetrial()'' needs.
-☛ Inspect the ''data_trl'' structure. Notice that there are now 88 trials instead of 1, and that each trial has its own time and data vectors. Plot the 47th trial as an example to verify this.
+☛ Inspect the ''data_trl'' structure. Notice that there are now 115 trials instead of 1, and that each trial has its own time and data vectors. Plot the 47th trial as an example to verify this.
 === Constructing the event-triggered spectrogram ===
@@ Line 357: / Line 327: @@
 cfg.foi          = 10:2:100; % frequencies of interest
 cfg.t_ftimwin    = ones(size(cfg.foi)).*0.5;  % window size: fixed at 0.5s
-cfg.toi          = -1:0.05:4; % times of interest
+cfg.toi          = -0.5:0.05:3.5; % times of interest
 TFR = ft_freqanalysis(cfg, data_trl);
@@ Line 370: / Line 340: @@
 You should get:
-{{ :analysis:course:week6_fig4.png?600 |}}
+{{ :analysis:nsb2014:week6_fig4.png?600 |}}
-This plot is an //average// spectrogram over 88 trials, with time 0 corresponding to the time of nosepoking into the reward receptacle. As with the examples using ''spectrogram()'', we can control the smoothness of the spectrogram by specifying the time and frequency steps. For the frequency steps this is accomplished the same way in FieldTrip, but for the time axis we can now simply specify a vector with the window centers instead of a window size and overlap.
+This plot is an //average// spectrogram over 115 trials, with time 0 corresponding to the time of nosepoking into the reward receptacle. As with the examples using ''spectrogram()'', we can control the smoothness of the spectrogram by specifying the time and frequency steps. For the frequency steps this is accomplished the same way in FieldTrip, but for the time axis we can now simply specify a vector with the window centers instead of a window size and overlap.
 Looking at the results, it appears there are some beta (~20Hz) and gamma events visible in the average spectrogram.
-☛ Change the ''cfg.toi'' and ''cfg.foi'' parameters to verify you can change the resolution of the time and frequency axes. What happens if you try to extend the time axis from -2 to 5?
+☛ Change the ''cfg.toi'' and ''cfg.foi'' parameters to verify you can change the resolution of the time and frequency axes. What happens if you try to extend the time axis from -2 to 5? What would you do if you wanted data for that extended range?
 === Options and parameters for FieldTrip spectrograms ===
@@ Line 405: / Line 375: @@
 You should now have a beautiful image that balances time and frequency resolution:
-{{ :analysis:course:week6_fig5.png?600 |}}
+{{ :analysis:nsb2014:ft_specgram2.png?600 |}}
 ☛ What causes the whitespace at the bottom of the image?
@@ Line 427: / Line 397: @@
 cfg.t_ftimwin    = 20./cfg.foi;  % 20 cycles per time window
-cfg.toi          = -2:0.05:0; % pre-nosepoke baseline
+cfg.toi          = -1:0.05:0; % pre-nosepoke baseline
 TFR_pre = ft_freqanalysis(cfg, data_trl);
-cfg.toi          = 0:0.05:2; % post-nosepoke
+cfg.toi          = 0:0.05:1; % post-nosepoke
 TFR_post = ft_freqanalysis(cfg, data_trl);
@@ Line 473: / Line 443: @@
 <code matlab>
+%% load the data
 fc = FindFiles('*.ncs'); % get filenames of all LFPs recorded
 data = ft_read_neuralynx_interp(fc); % load them all -- this will take a while
+data_all.hdr.Fs = data_all.fsample; % for some reason this is missing from the header
-% define layout
+%% define layout for later plotting
 cfg = [];
 cfg.layout = 'ordered'; cfg.channel = data.label;
@@ Line 487: / Line 459: @@
 <code matlab>
+%%
 cfg = [];
-cfg.trialfun = 'ft_trialfun_lineartracktone2';
+cfg.t = cat(2,getd(evt,'n0'),getd(evt,'n1'));
-cfg.trialdef.hdr = data.hdr;
+cfg.mode = 'nlx';
-cfg.trialdef.pre = 2.5;
+cfg.hdr = data_all.hdr;
-cfg.trialdef.post = 5;
+cfg.twin = [-1 4];
-cfg.trialdef.eventtype = 'nosepoke'; % could be 'nosepoke', 'reward', 'cue'
+trl = ft_maketrl(cfg);
-cfg.trialdef.location = 'both'; % could be 'left', 'right', 'both'
-cfg.trialdef.block = 'both'; % could be 'value', 'risk'
-cfg.trialdef.cue = {'c1','c3','c5'}; % cell array with choice of elements {'c1','c3','c5','lo','hi'}
-[trl, event] = ft_trialfun_lineartracktone2(cfg);
+cfg = [];
 cfg.trl = trl;
+data_trl = ft_redefinetrial(cfg,data_all);
-data_trl = ft_redefinetrial(cfg,data);
 %%
-cfg              = [];
+cfg              = []; % start with empty cfg
 cfg.output       = 'pow';
-%cfg.channel      = 'R016-2012-10-03-CSC04a';
 cfg.method       = 'mtmconvol';
 cfg.taper        = 'hanning';
-cfg.foi          = 1:1:100;
+cfg.foi          = 1:100; % frequencies of interest
-cfg.keeptrials = 'yes'; % should be needed for subsequent statistics...
+cfg.toi          = -0.5:0.05:3.5; % times of interest
-cfg.t_ftimwin    = 20./cfg.foi;  % 20 cycles per time window
+cfg.t_ftimwin = 20./cfg.foi;
-cfg.toi          = -1:0.05:4;
 TFR = ft_freqanalysis(cfg, data_trl);
 </code>
@@ Line 531: / Line 498: @@
 You should get:
-{{ :analysis:course:week6_fig6.png?600 |}}
+{{ :analysis:nsb2014:ft_specgram3.png?600 |}}
 This shows the baseline-corrected, event-aligned spectrograms for all 16 channels in this session. There is also an average shown (the subplot on the lower right). Note that a subset of the channels, presumably located in the ventral striatum, show a very similar time-frequency pattern, whereas another subset do not show this at all. The channels that do not show this same pattern are likely located in the hippocampus.
+==== Exercise ====
-==== Assignment ====
 Because the event-triggered spectrograms above show an average over trials, it is possible to lose touch with the properties of the raw data that generated it. For instance, the average may arise from a distribution of very similar looking individual trials, or it may be dominated by one atypical but extreme trial. For this reason it is important to compare the spectrogram to the raw data.
@@ Line 541: / Line 507: @@
 An example of a study where this was done well is [[http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3463873/ | Leventhal et al]]. Neuron 2012. [[http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3463873/figure/F3/ | Figure 3]] in this paper shows a number of event-triggered LFP traces, one for each trial.
-☛ Implement a function eventLFPplot(csc,event_times,varargin) to make such a plot, as follows:
+☛ Implement a function eventLFPplot(cfg,csc) to make such a plot. The ''cfg'' variable should contain an ''eventTimes'' field which is required, and optionally accept a ''window'' field to override the default window of [-1 3] seconds relative to the event times.
-<code matlab>
-% function eventLFPplot(csc,event_times,varargin)
-%
-% INPUTS
-%
-% csc: [1 x 1] mytsd, LFP signal to be plotted
-% event_times: [nEvents x 1] double with event times to align LFP on
-%
-% varargins (with defaults):
-%
-% t_window: [2 x 1] double indicating time window to use, e.g. [-1 3] for 1 second before to 3 seconds after event times
-</code>
 You can use a skeleton like the following:

User Tools

Site Tools

Differences

Page Tools