[FieldTrip] problem using ft_sourceparcellate

Fri Oct 8 17:05:31 CEST 2021

Hi Johannes,

On 8 Oct 2021, at 15:27, Johannes Gehrig <jgehrig at med.uni-frankfurt.de<mailto:jgehrig at med.uni-frankfurt.de>> wrote:

Dear Jan-Mathijs,

Thanks for your advice! It still took a while till Bob was my uncle, but after adding this (%JG) to ft_sourceparcellate in line 488 he was:

Well, I guess uncles in Germany are more often called Franz oder Peter...

function y = cellmean1(x)
siz = size(x);
if siz(1)==1 && siz(2)>1
  siz([2 1]) = siz([1 2]);
  x = reshape(x, siz);
end
if size(x,1)>0 %JG
    y = x{1};
    n = 1;
else %JG
    y = 1 %JG
    n = 1 %JG
end %JG

Using this I skipped the 5 parcells where there was no label at the "inside" compartment (->241 parcells).

It’s interesting that you need to make the above changes in order for the parcellation to work. I created the downsampled atlas was obtained as:

% load the standard 3d grid to-be-used
[ftver,path] = ft_version;
load(fullfile(path, 'template/sourcemodel', 'standard_sourcemodel3d8mm.mat'));
sourcemodel = ft_convert_units(sourcemodel, 'mm');

atlasfile = fullfile(path, 'template/atlas/brainnetome', 'BNA_MPM_thr25_1.25mm.nii');
atlas     = ft_read_atlas(atlasfile);

% downsample the atlas
cfg = [];
cfg.parameter = 'tissue';
cfg.interpmethod = 'nearest';
atlas_lowres = ft_sourceinterpolate(cfg, atlas, sourcemodel);
assert(numel(unique(atlas_lowres.tissue(:)))==numel(atlas.tissuelabel)+1);

and this seemed to work fine. That is (not shown here), a parcellation of a dummy dataset created on the fly worked without changing the cellmean function in ft_sourceparcellate.

This being said, you identified a new problem, which I didn’t think of before, and which you reference below:

Now the parcellation works, but the connectivity analyses is not possible.... I have the csd-, mom, and noisecsd- matrix but the connectivity analyses fails with:

"Error using getdimord (line 42)
field "crsspctrm" not present in data"

Do I have to label the mom matrix (241x3x1764 with dimord: chan_ori_rpttap) as crsspctrm and perform a svd analyses on the 2nd dimension?

Well indeed, the data bookkeeping in ft_connectivityanalysis (which causes the error for you) pre-dates the ft_sourceparcellate function, and I think we never got to making it more robust. Rather than fooling around with the dimensions, I’d like to suggest a different recipe (slightly different from the one you suggested below, although that also would have been a possibility). My suggested pipeline does not use ft_sourceparcellate, but rather uses ft_virtualchannel, in combination with 3 data objects: the lowres atlas, the sensor level freq-structure, and a source structure, obtained with ‘dics’, and keepfilter=‘yes’:

% alternative
cfg = [];
cfg.method = 'dics';
cfg.dics.keepfilter = 'yes';
% here you can add a  bunch of other options to be used, if needed
cfg.sourcemodel = leadfield;
source = ft_sourceanalysis(cfg, freq);

cfg = [];
cfg.parcellation = 'tissue';
cfg.numcomponent = 1;
freq_vc = ft_virtualchannel(cfg, freq, source, atlas_lowres);

cfg = [];
cfg.method = 'coh';
coh = ft_connectivityanalysis(cfg, freq_vc);

The leadfield variable contains precomputed leadfields on the same 8mm grid as the atlas. Now, one thing that one can play with is the numcomponent option for ft_virtualchannel. The long story short is that the parcellation step by default involves a SVD (~PCA) decomposition of the concatenated dipoles’ data for all dipoles that belong to a single parcel. This may be increased at libitum, but this will yield a connectivity matrix that does not directly map one-to-one on the individual parcel pairs anymore. For now, I think I would recommend using numcomponent = 1;

I was thinking if not maybe tlock data -> lcmv beamformer -> parcellation -> freq.analysis -> connectivity analysis is the "better" option?

Yes, this would have been a possibility in principle, but tricky, since you’d need single trial data at the source level for the spectral analysis step, so Franz or Peter would have to wait a bit on this route, too.

Best wishes,
JM

Thanks one again for your support and best regards,
Johannes

Am 07/10/2021 um 14:45 schrieb Schoffelen, J.M. (Jan Mathijs):
Dear Johannes,

Gotcha, thanks for the information. I think that in general it would be wise to use a volumetric atlas since you are using a volumetric grid for source reconstruction. For sure a 3D version of the Glasser atlas exists, and it should not be too difficult to trick FieldTrip into importing it. Do you have a specific reason for using this particular parcellation, or would any other sensible number-of-points reducing atlas do equally well for you? If the latter is true, you could use the brainnetome atlas (which is shipped with fieldtrip, so easy to use).

Then, under the assumption that you are using inverse warped 3D grid for the patients, it would make most sense to create a low-resolution version of the atlas, sampled on the standard grid that you use for the beamformer analysis: something like ft_sourceinterpolate(cfg, atlas, sourcemodel), where the atlas is the atlas (in normalised MNI space), the sourcemodel is the 8mm standard sourcemodel, and cfg.method = ’nearest’.

Ideally in the low-res version of the atlas each dipole position in the ‘inside’ compartment should be labeled (and vice versa each label in the original high res atlas should be present in the low res atlas). If these constraints are not met, the parcellation might still fail later on.

Then, once you are ready for the parcellation, you need to ensure the source.pos and atlas (low-res).pos to be identical, and Bob should be your uncle.

There might be some dedicated documentation about this on the website, but since it is currently down, I couldn’t search for and point you to it.

Best wishes,
Jan-Mathijs

On 7 Oct 2021, at 12:02, Johannes Gehrig <jgehrig at med.uni-frankfurt.de<mailto:jgehrig at med.uni-frankfurt.de>> wrote:

Dear Jan-Mathijs,

In general, I want to compare the connectivity between two MEG sessions in which patients took different drugs.

I followed the tutorial "whole brain connectivity" to get a first idea of the data (so I tried to calculate the connectivity in dipole space...).
I tried to use the parcellation of Glasser et al 2016 (atlas_MMP1.0_4k.mat). To fit this to my data, I used ft_sourceinterpolate of my source model (subject-specific, using the template "standard_sourcemodel3d8mm") and the atlas.

I think this is exactly the cause: "One possible cause is the fact that some parcels in the input parcellation are not assigned any of the locations where your subject data is represented."
Because when I debug (ft_sourceparcellate), the second parcel is not found (j==2, and dat(tissue==j) is an "Empty cell array: 0-by-1").

How could I solve this issue?

Thanks in advance an best regards,
Johannes

Von: Schoffelen,J.M. (Jan Mathijs) via fieldtrip<mailto:fieldtrip at science.ru.nl>
Gesendet: Donnerstag, 7. Oktober 2021 10:28
An: FieldTrip discussion list<mailto:fieldtrip at science.ru.nl>
Cc: Schoffelen, J.M. (Jan Mathijs)<mailto:janmathijs.schoffelen at donders.ru.nl>
Betreff: Re: [FieldTrip] problem using ft_sourceparcellate

Dear Johannes,

Could you explain a bit more what you attempt to do? Specifically: what is the parcellation that you are using? Ideally, you should use one that is defined on exactly the same surface as your data.  Otherwise (e.g. when using a parcellation defined in 3D space) some on-the-fly interpolation is required, and this may not work robustly enough. One possible cause is the fact that some parcels in the input parcellation will not get assigned any of the locations at which your functional data is represented. Which is probably the cause of the error you are facing.

Best wishes,
Jan-Mathijs

PS: as a side note, in general I would recommend to compute connectivity in parcellated space, rather than in ‘dipole space’. Apart from potential memory issues, the sheer number of pairwise combination makes the interpretation / statistical evaluation quite complicated.

On 7 Oct 2021, at 09:08, Johannes Gehrig via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>> wrote:

Dear Stefan,

thanks for your reply. At the moment I use simple "coherence" following the „whole brain connectivity“ tutorial. However, if there was a memory efficient solution for dwPLI that is also an option. If I followed the discussion correctly, this is not yet fully implemented?

Have you tried to parcellate your data before doing the connectivity analysis? I suspect that my problem could be related with ft_sourceinterpollate, because the output of the source analysis (avg.csd, avg.mom, avg.noisecsd, computed withels pcc) is no longer present after ft_sourceinterpolate and without the source interpolation the source parcellation doesn't work, so far my impression...

Thanks and best regards,
Johannes

Von: Stefan Dvoretskii<mailto:stefan.dvoretskii at tum.de>
Gesendet: Mittwoch, 6. Oktober 2021 17:39
An: FieldTrip discussion list<mailto:fieldtrip at science.ru.nl>
Cc: Johannes Gehrig<mailto:jgehrig at med.uni-frankfurt.de>; Gil Avila, Cristina<mailto:cristina.gil at tum.de>
Betreff: Re: [FieldTrip] problem using ft_sourceparcellate

Dear Johannes,

which connectivity measure are you using? There is a potential bottleneck in the 'ft_connectivityanalysis', which can however be circumvented with custom code, at least for dwPLI. Then you won't even need to parcelate.
We are currently working to incorporate it into the Fieldtrip code.

Best regards,
Stefan

пн, 4 окт. 2021 г., 11:06 Johannes Gehrig via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>>:
Dear Fieldtrip-user,
I am facing a problem using ft_sourceinterpolate. I am following the "Whole brain connectivity and network analysis"- Tutorial and as I can not compute the connectivity directly with the source data ("out of memory") I would like to perform the source parcellation before computing the connectivity. Unfortunately (fieldtrip 20211001 and matlab2018b) ft_sourceparcellate crashes after calling it. Before computing the source data the sourcemodel was interpolated with the atlas using ft_sourceinterpolate. Therefore the amount of positions in the atlas and the source data is the same (8004).
load('atlas_MMP1.0_4k.mat');
cfg = [];
cfg.parcellation = 'parcellation';
cfg.method       = 'mean';
cfg.parameter    = 'all';
parc_sourceOn = ft_sourceparcellate(cfg, source_ON, atlas);
Error message:"
there are in total 8004 positions, 4544 positions are inside the brain, 8004 positions have a label
4544 of the positions inside the brain have a label
4544 of the labeled positions are inside the brain
0 of the positions inside the brain do not have a label
creating 362 parcels for parameter csd by taking the mean
computing parcellation
232           tmp(j,:,:) = cellmean1(dat(tissue==j));
computing parcellation for L_V1_ROI
Index exceeds array bounds.

Error in ft_sourceparcellate>cellmean1 (line 494)
y = x{1};

Error in ft_sourceparcellate (line 232)
          tmp(j,:,:) = cellmean1(dat(tissue==j));
"
Any help and ideas are very much appreciated!
Best regards, Johannes
_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://doi.org/10.1371/journal.pcbi.1002202
<https://urldefense.com/v3/__https:/doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!pg0W3nFXzJCE9mE1pvYfYPxllCL06pTokAzKEH8cYGsjnXbtXGbXLrC_6SLHcu0eFLPB1lmxevpST48$>

_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://urldefense.com/v3/__https://doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!pg0W3nFXzJCE9mE1pvYfYPxllCL06pTokAzKEH8cYGsjnXbtXGbXLrC_6SLHcu0eFLPB1lmxevpST48$<https://urldefense.com/v3/__https:/doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!pg0W3nFXzJCE9mE1pvYfYPxllCL06pTokAzKEH8cYGsjnXbtXGbXLrC_6SLHcu0eFLPB1lmxevpST48$>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20211008/36725387/attachment.htm>