Hi again Giles,
by mistake I sent you an intermediate version of the m-files.
Please ignore the scripts from the previous e-mail and please find the
proper version attached here.
Best
Giorgos
On 12/02/2015 12:48, Giles Colclough wrote:
Hi,
I have two queries I'm looking for help with.
1. Why is there a scaling difference between the magnitude of the data
recorded in the noise scans, and the output of the rmegpreproc pipeline?
I find about 30 orders of magnitude difference between the variance of
the empty room scan and the variance of the data outputted by
rmegpreproc.
Have the data been uniformly scaled? If so, by how much? This would be
useful information, for example, when source localising using the
empty room data as an estimate of the noise.
2. When source localizing with a beamformer, I find an unusual
variance profile in the resting state data.
Normally when we look at resting data, the source variance is
concentrated in a 'halo' around the outside of the cortex. In the HCP
data, this halo is present, but there is also a bright stripe
bisecting the brain around the central sulcus. We find this suprising,
but have been unable to determine what's causing it.
The effect is very repeatable between participants.
I attach a screenshot of the effect, and two scripts to replicate the
analysis.
Any help or insights would be greatly appreciated.
Many thanks,
Giles
_______________________________________________
HCP-Users mailing list
HCP-Users@humanconnectome.org
http://lists.humanconnectome.org/mailman/listinfo/hcp-users
---
This email is free from viruses and malware because avast! Antivirus protection
is active.
http://www.avast.com
_______________________________________________
HCP-Users mailing list
HCP-Users@humanconnectome.org
http://lists.humanconnectome.org/mailman/listinfo/hcp-users
%% PART 1 - SENSOR LEVEL - COMPARE POWER and VARIANCE BETWEEN RESTING STATE AND
ROOM NOISE
%---------------------
% GIORGOS: Add fieldtrip to the path
%{
fieldtripdir='/home/michalareasg/toolboxes/matlab/fieldtrip/fieldtrip_svn';
addpath(fieldtripdir);
ft_defaults;
%}
%% Load in relevant data
datadir = '/mnt/hps/slurm/michalareasg/pipeline/TestGilEmail/';
subjID = 185442;
iSession = 5;
megDir = fullfile(datadir, num2str(subjID), '/MEG/Restin/rmegpreproc/');
anatomyDir = fullfile(datadir, num2str(subjID), '/MEG/anatomy/');
noiseDir = fullfile(datadir, num2str(subjID),
'/unprocessed/MEG/1-Rnoise/4D/');
megFile = fullfile(megDir, sprintf('%d_MEG_%d-Restin_rmegpreproc.mat',
subjID, iSession));
noiseFile = fullfile(noiseDir,'c,rfDC');
headFile = fullfile(anatomyDir, sprintf('%d_MEG_anatomy_headmodel.mat',
subjID));
sourceFile = fullfile(anatomyDir,
sprintf('%d_MEG_anatomy_sourcemodel_3d8mm.mat', subjID)); % use 8mm resolution
%Load Raw Noise data
%===============================================================
% Split Room Noise in pseudo trials of 2 seconds
noiseHdr=ft_read_header(noiseFile);
pseudoDur=2; %pseudo Trial Duration
pseudoDurSamps=ceil(pseudoDur*noiseHdr.Fs);
pseudoTrl=[1:pseudoDurSamps:noiseHdr.nSamples-pseudoDurSamps]';
pseudoTrl=[pseudoTrl pseudoTrl+(pseudoDurSamps-1)];
pseudoTrl=[pseudoTrl zeros(size(pseudoTrl,1),1)];
cfg = [];
cfg.channel='MEG';
cfg.dataset = noiseFile;
cfg.trl=pseudoTrl;
orignoisedata=ft_preprocessing(cfg);
%------------------------------------------------------------------------------------------
% Detrend each pseudo trial just in case there is any problematic channel with
any
% significant slow increasing or decreasing trends. Should make little
% difference for channels with no problems
cfg=[];
cfg.detrend='yes';
noisedata=ft_preprocessing(cfg,orignoisedata);
%============================================================================
% LOAD RESTING DATA
load(megFile); % loads resting data data
%============================================================================
% select the common channels between resting and noise
[indA,indB]=match_str( data.label, noisedata.label);
noisedata=ft_selectdata(noisedata,'channel',indB);
%============================================================================
% CONCATANATE TRIALS
datRest=cat(2, data.trial{:});
datNoise=cat(2, noisedata.trial{:});
%============================================================================
% COMPUTE COVARIANCE
covRest=cov(datRest');
covNoise=cov(datNoise');
varRest=std(datRest,0,2).^2;
varNoise=std(datNoise,0,2).^2;
% PLOT VARIANCE (the diagonal of the Covariance Matrix)
figure; hold on
plot(diag(covNoise),'r');
plot(diag(covRest),'g');
ylabel('power');
xlabel('channel index');
legend('Room Noise',' Resting State','location','northwest');
%===============================================================
% GIORGOS: In the above plot you will see that the variance of all channels is
very
% similar apart from one channel (index=211, label='A197') which
% appears to have many jumps during the noise measurement but is stable
% during the resting state.
% You can check this channels abnormal behavior during noise scan by
% ft_databrowser([],orignoisedata);
% But even for this channel the variance are only about one order
% of magnitude higher than the rest of the channels
%% PART 2 - SOURCE LEVEL - THE EFFECT OF NORMALIZATION
%% LOAD SOURCE SPACE
load(headFile); % loads headmodel
load(sourceFile); % loads sourcemodel3d
%% Get Gradiometer definition from data
grad = data.grad;
% GILES:
% apply gradiometers - copied from hcp_icamne L100.
% grad = ft_apply_montage(grad, grad.balance.Supine, 'keepunused', 'yes',
'inverse', 'yes');
% GIORGOS: YOU DONT NEED TO DO THE ABOVE BECAUSE THE GRADIOMETER HAS ALREADY THE
% SUPINE BALANCING REMOVED.
% Your grad should have the following fields
% grad.balance=
% Supine: [1x1 struct]
% pca: [1x1 struct]
% current: 'invcomp'
% previous: {'pca' 'none'}
% invcomp: [1x1 struct]
% meaning that first the principal components('pca') of the reference
% sensors which measure enviromental noise were removed from data and
% then the Independent Components from heart and eye
% artifacts('invcomp') were removed. So in the balancing there are not
% SUppine balnce weights involved so you dont need to do the inverse
% Supine balancing.
%% Convert everything into mm
sourcemodel3d = ft_convert_units(sourcemodel3d,'mm');
headmodel = ft_convert_units(headmodel,'mm');
grad = ft_convert_units(grad,'mm');
%% Compute the leadfields
cfg = [];
cfg.vol = headmodel;
cfg.grid = sourcemodel3d;
cfg.grad = grad;
cfg.channel = data.label;
cfg.normalize = 'no';
cfg.reducerank = 2;
gridLF = ft_prepare_leadfield(cfg);
% GIORGOS: Above you select cfg.reducerank = 2 so the resulting leadfield
% has already a reduced rank.
%% Reshape data
% sourcedata = data; % GIORGOS : Not used anywhere so commented out
dat = cat(2, data.trial{:});
time = (1:size(dat,2))./data.fsample; % not real time, due to bad
segments having been removed
Fs = data.fsample;
%% Bandpass filter before source localization
% use wide-band 4-30 Hz reconstruction
Fbp = [4 30];
N = 5;
filteredData = ft_preproc_bandpassfilter(dat, data.fsample, Fbp, N);
%% Compute the beamformer weights
C = cov(dat.');
% reduce dimensionality
dimensionality = 220;
invC = pinv_plus(C,dimensionality);
% rank for leadfields
reduce_rank = 2;
%% GIORGOS: CASE A: ORIGINAL GILES CODE - BEAMFORMER SPATIAL WEIGHTS ARE
NORMALISED BY THEIR NORM
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% Apply weights normalisation
W{i} = W{i} ./ sqrt(W{i}*W{i}');
% GIORGOS: I am not sure why here you normalise the Spatial Filter with
% its norm. This makes the product W{i}*lf ~= 1 as it should be
% according to beamformer defintion.
%
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
f1= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - Giles code - Beamf. Weights normalised by their Norm');
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score to
avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
ft_sourceplot(cfg,sourceloc);
%% GIORGOS: CASE B: NO NORMALIZATION IN BEAMFORMER SPATIAL WEIGHTS , NOR IN
THE LEADFIELD
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% GIORGOS: As here there is not beamformer weight normalization , but instead
the beamformer weights are computed based on the normalizzed leadfield, the
equation W{i}*lf = 1 hold true as it should.
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
f2= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - No normalization of Beamf. Weights NOR of Leadfield');
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score to
avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
ft_sourceplot(cfg,sourceloc);
%% GIORGOS: CASE C: NORMALIZATION OF LEADFIELDS
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
%----------------------------------------
lf_mag(i) = norm(lf);
lf= lf ./ lf_mag(i);
% The leadfields have an inherent bias towards the center of the brain.
% This because a source at the center of the brain needs much more
% power to produce the same signal at the sensors as compared to a
% source located on the cortex. The typical way to remove this bias is
% to normalise the LEADFIELDS by their norm (and not the Beamformer
% Spatial Weights);
%----------------------------------------
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% GIORGOS: As here there is not beamformer weight normalization , but
instead the beamformer weights are computed based on the normalizzed leadfield,
the equation W{i}*lf = 1 hold true as it should.
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
f3= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - Leadfield Normalization with its Norm');
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score to
avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
ft_sourceplot(cfg,sourceloc);
%% FINAL COMMENT :
%%========================================================================
% The leadfields have an inherent bias towards the center of the brain.
% This because a source at the center of the brain needs much more
% power to produce the same signal at the sensors as compared to a
% source located on the cortex. The typical way to remove this bias is
% to normalise the LEADFIELDS by their norm (and not the Beamformer
% Spatial Weights);
% so for the normalisation lf is basically divided by
% (l(1).^2+l(2).^2+l(3).^2+.......+l(Nsensors).^2)^(0.5) which is the
% norm.
% the power exponent 0.5 corresponds to the square root.
% If this exponent is 0 , corresponds to no normalization. The higher the
% exponent the more the bias is reversed from the centre to the surface.
%% GIORGOS: CASE D: NORMALIZATION OF LEADFIELDS
for normalizeparam = [0:0.25:1.5],
normalizeparam
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
%----------------------------------------
lf_mag(i) = sum(lf(:).^2)^normalizeparam;
lf= lf ./ lf_mag(i);
% The leadfields have an inherent bias towards the center of the brain.
% This because a source at the center of the brain needs much more
% power to produce the same signal at the sensors as compared to a
% source located on the cortex. The typical way to remove this bias is
% to normalise the LEADFIELDS by their norm (and not the Beamformer
% Spatial Weights);
%----------------------------------------
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% GIORGOS: As here there is not beamformer weight normalization , but
instead the beamformer weights are computed based on the normalizzed leadfield,
the equation W{i}*lf = 1 hold true as it should.
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
%{
f3= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - Leadfield Normalization with its Norm');
%}
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score
to avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
cfg.colormap='hot';
ft_sourceplot(cfg,sourceloc);
end
%% GIORGOS COMMENT:
% From the figures plotted here you can see the bias moving from the center to
the periphery
% For the exponent value of 0.5 the power seems similar to what you showed
% in your plot. THe normalization by the norm (i.e. 0.5) creates an overall bias
% profile, in intermediate depth between the centre and the surface. In
% order to bring the maxima closer to the surface , a exponent of 1 is more
% appropriate
