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
% script to replicate weird source variance activity in HCP MEG data
% Giles Colclough
% 12 February 2015
%% Load in relevant data
dataDir = '/home/gileslc/data/hcp/HCPsensor/';
subjID = 185442;
iSession = 5;
megDir = fullfile(datadir, num2str(subjID), '/MEG/Restin/rmegpreproc/');
anatomyDir = fullfile(datadir, num2str(subjID), '/MEG/anatomy/');
megFile = fullfile(megDir, sprintf('%d_MEG_%d-Restin_rmegpreproc.mat', subjID, iSession));
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(megFile); % loads data
load(headFile); % loads headmodel
load(sourceFile); % loads sourcemodel3d
mni_brain_2mm = '/home/gileslc/OSL-Repo/osl2/std_masks/MNI152_T1_2mm_brain.nii.gz';
mni_brain_8mm = '/home/gileslc/OSL-Repo/osl2/std_masks/MNI152_T1_8mm_brain.nii.gz'; % downsampled MNI152_T1_2mm_brain.nii.gz
%% Prepare the data
grad = data.grad;
% apply gradiometers - copied from hcp_icamne L100.
grad = ft_apply_montage(grad, grad.balance.Supine, 'keepunused', 'yes', 'inverse', 'yes');
% 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;
% prepare leadfields
gridLF = ft_prepare_leadfield(cfg);
%% Reshape data
sourcedata = data;
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;
for i = length(gridLF.inside):-1:1;
lf = gridLF.leadfield{gridLF.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);
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}');
end%for
% reformat W
weights = cat(1, W{:});
%% Compute source variance
sourceVar = diag(weights * C * weights.');
%% Output nifti file
% we need to to a bit of convoluted resampling to get into mni space for
% fslview
% declare empty volume of same size as grid
vol = zeros(length(sourcemodel3d.cfg.grid.template.xgrid), ...
length(sourcemodel3d.cfg.grid.template.ygrid), ...
length(sourcemodel3d.cfg.grid.template.zgrid));
% fill in volume with source variance at correct locations
for vox = sourcemodel3d.cfg.grid.template.inside;
vol(sourcemodel3d.cfg.grid.template.pos(vox,1) == sourcemodel3d.cfg.grid.template.xgrid,...
sourcemodel3d.cfg.grid.template.pos(vox,2) == sourcemodel3d.cfg.grid.template.ygrid,...
sourcemodel3d.cfg.grid.template.pos(vox,3) == sourcemodel3d.cfg.grid.template.zgrid) = ...
sourceVar(vox == sourcemodel3d.cfg.grid.template.inside);
end%for
% resample to appropriate nifti space
[~, dims, scales] = read_avw(mni_brain_8mm);
mni_grid_x = 90 - (0:dims(1)-1)*scales(1);
mni_grid_y = -126 + (0:dims(2)-1)*scales(2);
mni_grid_z = -72 + (0:dims(3)-1)*scales(3);
% turn grid vectors into full matrices
[mni_grid_x, mni_grid_y, mni_grid_z] = ndgrid(mni_grid_x, mni_grid_y, mni_grid_z);
x_mm = round(sourcemodel3d.cfg.grid.template.xgrid * 10);
y_mm = round(sourcemodel3d.cfg.grid.template.ygrid * 10);
z_mm = round(sourcemodel3d.cfg.grid.template.zgrid * 10);
% perform interpolation
Vq = interpn(x_mm, y_mm, z_mm,...
vol, ...
mni_grid_x, ...
mni_grid_y, ...
mni_grid_z);
% output nifti
niiFileName = 'source_variance.nii';
niiFileName_2mm = 'source_variance_ds2mm.nii';
save_avw(Vq, niiFileName, 'f', [8 8 8 1]);
% dilate volume to avoid edge effects
unix(['fslmaths ' niiFileName ' -kernel gauss 8 -dilM ' niiFileName_2mm]);
% resample
unix(['flirt -in ' niiFileName_2mm ' -applyxfm -init ' getenv('FSLDIR') '/etc/flirtsch/ident.mat -out ' niiFileName_2mm ' -paddingsize 0.0 -interp trilinear -ref ' mni_brain_2mm]);
% mask
unix(['fslmaths ' niiFileName_2mm ' -mas ' mni_brain_2mm ' ' niiFileName_2mm]);
%% View
% look for an unusual high variance plane bisecting the front and back of
% the brain - about the location of the central sulcus
% This is markedly different to the normal 'halo' of high power around the
% outside of the brain
unix(['fslview ' mni_brain_8mm ' ' niiFileName ' -l Red-Yellow -b 6e-27,1.5e-26 -t 0.9 &'])
unix(['fslview ' mni_brain_2mm ' ' niiFileName_2mm ' -l Red-Yellow -b 6e-27,1.5e-26 -t 0.9 &'])
function [X, pca_order] = pinv_plus(A,varargin)
% PINV_PLUS Modified Pseudoinverse which allows a specified
% dimensionality (order) to be passed in.
%
% pinv_plus(A) is the same as pinv(A)
%
% [X, pca_order]=pinv_plus(A,pca_order) uses at most only the top pca_order PCs to compute inverse
%
% pinv_plus(A,pca_order,1) uses precisely the top pca_order PCs to compute inverse
%
% If pca_order=-1 then calls spm_pca_order to compute order to use
%
% Class support for input A:
% float: double, single
%
% Mark Woolrich
tol=[];
r=size(A,1);
if nargin >= 2,
pca_order = varargin{1};
if(pca_order==-1),
pca_order=spm_pca_order(A);
end;
else
pca_order = size(A,1);
end;
max_r=pca_order;
if nargin >= 3,
fix_to_max_r=varargin{2};
else
fix_to_max_r=0;
end;
if isempty(A) % quick return
X = zeros(size(A'),class(A));
return
end
% [U,S,V] = svd(A,0); s = diag(S);
[m,n] = size(A);
if n > m
X = pinv(A',varargin{:})';
else
if(max_r<size(A,1))
%[U,S,V] = svds(A,max_r);
[U,S,V] = svd(A,0);
max_r=min(length(diag(S)),max_r);
else
[U,S,V] = svd(A,0);
end;
if m > 1, s = diag(S);
elseif m == 1, s = S(1);
else s = 0;
end
tol=[];
if max_r < 0
diffs=abs(diff(s)./s(1:end-1));
for i=2:length(diffs),
if(diffs(i) > 6*mean(diffs(1:i))),
break;
end;
end;
max_r = i;
%figure;plot(diffs);ho;plot(r,diffs(r),'*');
end;
% max_r
% fix_to_max_r
if(fix_to_max_r==1),
r=max_r;
else,
tol=max(m,n) * eps(max(s));
%tol=max(size(A)) * norm(A) * eps(class(A));
r = sum(s > tol);
r = min(r,max_r);
end;
%figure;plot(s);
s2=s;
if (r == 0)
X = zeros(size(A'),class(A));
else
s = diag(ones(r,1)./s(1:r));
X = V(:,1:r)*s*U(:,1:r)';
end
pca_order=r;
end
