Sourcespace Analysis
This tutorial guides you through the process of performing sourcespace analysis using Beamforming in Fieldtrip. For more background on beamformers, please see this introduction by Hilliebrand & Barnes (2005).
We will be analysing the same data as analysed in the preprocessing tutorial, involving the presentation of a high contrast visual grating for 1.5s.
Set-up Paths
The first step will be set up the paths on your machine.
You will need to following toolboxes, which can be downloaded from Github
% Path to the Fieldtrip Toolbox
path_to_fieldtrip = '/Users/rseymoue/Documents/GitHub/fieldtrip/';
% Path to MQ_MEG_Scripts
path_to_MQ_MEG_Scripts = '/Users/rseymoue/Documents/GitHub/MQ_MEG_Scripts/';
% Add Fieldtrip
addpath(path_to_fieldtrip);
ft_defaults;
% Add MQ_MEG_Scripts
addpath(genpath(path_to_MQ_MEG_Scripts));
Load Variables Required for Beamforming
We need to load the preprocessed MEG data (data_clean.mat) and singleshell headmodel, sourcemodel and transformed MEG sensors made during MRI-MEG coregistration and by mq_realign_sens.
% Load preprocessed MEG data
load('data_clean.mat');
% Load singleshell headmodel
load('headmodel.mat');
% Load sourcemodel
load('sourcemodel3d.mat');
% Load transformed MEG sensors
load('grad_trans.mat');
Band-pass Filter and Select Time-points of Interest
Here we are interested in localising the increase in gamma-band (40-70Hz) power for time-points between 0.3 to 1.5s post-stimulus onset; versus a 1.2s baseline period from -1.5s to -1.3s.
%% BP Filter
cfg = [];
cfg.bpfilter = 'yes';
cfg.bpfreq = [40 70]; %band-pass filter in the required range
data_filtered = ft_preprocessing(cfg,data_clean);
%% Here we redefine trials based on the time-points of interest.
% Make sure the timepoints are of equivalent length
cfg = [];
cfg.toilim = [-1.5 -0.3];
datapre = ft_redefinetrial(cfg, data_filtered);
cfg.toilim = [0.3 1.5];
datapost = ft_redefinetrial(cfg, data_filtered);
Prepare Leadfield
By using ft_prepare_leadfield, we can compute the forward model for many dipole locations on a 3D sourcemodel and store it for efficient inverse modelling.
%% Create leadfield
cfg = [];
cfg.channel = data_clean.label;
cfg.sourcemodel = sourcemodel3d;
cfg.headmodel = headmodel;
cfg.grad = grad_trans;
cfg.normalize = 'yes';
lf = ft_prepare_leadfield(cfg);
% make a figure of the single subject{i} headmodel, and grid positions
figure; hold on;
ft_plot_vol(headmodel, 'facecolor', 'cortex', 'edgecolor', 'none');
alpha 0.5; camlight;
ft_plot_mesh(lf.pos(lf.inside,:));
ft_plot_sens(grad_trans, 'style', 'r*')
Compute Covariance Matrix
Fairly self-explanatory!
We also create an AVERAGE covariance matrix for the common spatial filter
%% Compute Covariance Matrix
% ...using ft_timelockanalysis
cfg = [];
cfg.channel = data_clean.label;
cfg.covariance = 'yes';
avgpre = ft_timelockanalysis(cfg,datapre);
avgpst = ft_timelockanalysis(cfg,datapost);
% Here we are making an average covariance matrix from the pre and post
% time periods to make a common filter
avg = avgpst;
avg.cov = (avgpst.cov + avgpre.cov)/2;
Source Analysis
Now we can call ft_sourceanalysis with the follwing cfg options:
%% Do ya beamforming
cfg = [];
cfg.channel = data_clean.label;
cfg.method = 'lcmv';
cfg.sourcemodel = lf;
cfg.headmodel = headmodel;
cfg.lcmv.keepfilter = 'yes';
cfg.lcmv.lambda = '5%';
cfg.lcmv.fixedori = 'yes';
cfg.grad = grad_trans;
sourceavg = ft_sourceanalysis(cfg, avg);
% Now do beamforming for the two time points separately using the same spatial
% filter computed from the whole trial
cfg.sourcemodel.filter = sourceavg.avg.filter; %uses the grid from the whole trial average
%Pre-grating
sourcepreS1 = ft_sourceanalysis(cfg, avgpre);
%Post-grating
sourcepstS1 = ft_sourceanalysis(cfg, avgpst);
Make sure your field positions match the template grid
In order for the source positions to be consistent across participants, we need to make the .pos field equal to the template (MNI) grid. For more information on this, plese see: LINK
% Load template sourcemodel (8mm)
load([path_to_fieldtrip 'template/sourcemodel/'...
'standard_sourcemodel3d8mm.mat']);
template_grid = sourcemodel;
template_grid = ft_convert_units(template_grid,'mm');
clear sourcemodel;
sourcepreS1.pos = template_grid.pos;
sourcepstS1.pos = template_grid.pos;
Calculate Percentage Change
Here we use ft_math to keep track of things.
cfg = [];
cfg.parameter = 'avg.pow';
cfg.operation = '((x1-x2)./x2)*100';
sourceR = ft_math(cfg,sourcepstS1,sourcepreS1);
Interpolate onto SPM T1 Standard Brain
Interpolate back to the SPM stnadard T1 brain
% Load the SPM T1 brain
mri = ft_read_mri([path_to_fieldtrip 'template/anatomy/'...
'single_subj_T1.nii']);
cfg = [];
cfg.voxelcoord = 'no';
cfg.parameter = 'pow';
cfg.interpmethod = 'nearest';
sourceI = ft_sourceinterpolate(cfg, sourceR, mri);
Produce a 2D Plot
Try playing around with the settings of ft_sourceplot
%% Produce 2D Plot
cfg = [];
cfg.funparameter = 'pow';
ft_sourceplot(cfg,sourceI);
ft_hastoolbox('brewermap', 1); % ensure this toolbox is on the path
colormap(flipud(brewermap(64,'RdBu'))) % change the colormap
Produce a 3D Plot
Try playing around with come of the cfg options!
%% Produce 3D Plot
cfg = [];
cfg.method = 'surface';
cfg.funparameter = 'pow';
cfg.projmethod = 'nearest';
cfg.surfinflated = 'surface_inflated_both.mat';
%cfg.surfdownsample = 10
cfg.projthresh = 0.2;
cfg.camlight = 'no';
ft_sourceplot(cfg, sourceI);
ft_hastoolbox('brewermap', 1); % ensure this toolbox is on the path
colormap(flipud(brewermap(64,'RdBu'))) % change the colormap
light ('Position',[0 0 50])
light ('Position',[0 -50 0])
material dull;
drawnow;
view([0 0]);
title('% Change Gamma Power vs Baseline');
set(gca,'FontSize',14);
Create Two Virtual Electrodes Based on the AAL Atlas
Here is some code to create two 'virtual electrodes' or 'virtual channels'. We use the AAL atlas, with ROIs in the left calcarine and right calcarine sulcus.
More technically, this code concatenates all spatial filters within a ROI, multipled by sensor-level data, followed by an SVD. By extracting the first principle component we end up with just ONE filter for each ROI, which we can right-multiply with the sensor-level data.
This is followed by some time-frequency analysis.
% Read in the AAL atlas from Fieldtrip and
atlas = ft_read_atlas([path_to_fieldtrip 'template/atlas/aal/ROI_MNI_V4.nii']);
atlas = ft_convert_units(atlas,'mm');
disp('Creating ROIs in L_Calcarine and R_Calcarine');
% Interpolate the atlas onto 8mm grid
cfg = [];
cfg.interpmethod = 'nearest';
cfg.parameter = 'tissue';
sourcemodel2 = ft_sourceinterpolate(cfg, atlas, template_grid);
labels = {'Calcarine_L','Calcarine_R'};
for lab = 1:length(labels)
% Find atlas points
atlas_points = find(sourcemodel2.tissue==...
find(contains(sourcemodel2.tissuelabel,labels{lab})));
% Concatenate spatial filter
F = cat(1,sourceavg.avg.filter{atlas_points});
% Do SVD
[u,s,v] = svd(F*avg.cov*F');
filter = u'*F;
% Create VE
VE = [];
VE.label = {labels{lab}};
VE.trialinfo = data_clean.trialinfo;
for subs = 1:(length(data_clean.trialinfo))
% note that this is the non-filtered "raw" data
VE.time{subs} = data_clean.time{subs};
VE.trial{subs}(1,:) = filter(1,:)*data_clean.trial{subs}(:,:);
end
% Save to file with label and subject information
save(sprintf('VE_%s',labels{lab}),'VE');
% Compute TFR using multitapers
cfg = [];
cfg.method = 'mtmconvol';
cfg.output = 'pow';
cfg.pad = 'nextpow2'
cfg.foi = 1:1:100;
cfg.toi = -2.0:0.02:3.0;
cfg.t_ftimwin = ones(length(cfg.foi),1).*0.5;
cfg.tapsmofrq = ones(length(cfg.foi),1).*8;
multitaper = ft_freqanalysis(cfg, VE);
% Plot
cfg = [];
cfg.ylim = [30 100];
cfg.xlim = [-0.5 1.5];
cfg.baseline = [-1.5 -0.3];
cfg.baselinetype = 'db';
cfg.zlim = 'maxabs';
figure; ft_singleplotTFR(cfg, multitaper);
title(labels{lab},'Interpreter','none');
ft_hastoolbox('brewermap', 1); % ensure this toolbox is on the path
colormap(flipud(brewermap(64,'RdBu'))) % change the colormap
xlabel('Time (s)'); ylabel('Frequency (Hz)');
set(gca,'FontSize',20); drawnow;
print(sprintf('VE_TFR_%s',labels{lab}),'-dpng','-r200');
end
Group Statistical Analysis
TO FOLLOW...