[FieldTrip] ft_sourceanalysis single trial reconstruction
Schoffelen, J.M. (Jan Mathijs)
janmathijs.schoffelen at donders.ru.nl
Sat Oct 30 19:07:23 CEST 2021
Dear Manon,
Just as a pre-preliiminary: you write that you want to compute EMG-MEG coherence. Did you consider to use DICS for this? If you know which frequency (range) you want to look at, I think that this would be the most straightforward. There’s plenty of documentation on the website for this.
Then, as a second preliminary: you write that you want to use a ‘beamformer technique in the time domain’, but in your code you use ‘mne’ as a method, which is a different algorithm altogether. Yet, it is not really relevant to the problem you are facing, so for now I’ll let it slip :). As a side note, the option ‘rawtrial’, for ft_sourceanalysis, and as far as I know, does not have any funcionality for the mne method, and only (historically) was functional for beamformers. These days, this option is explicitly DISrecommended. (nobody has found the time yet to cleanly deprecate this functionality altogether)
OK, after these preliminaries, let’s think about the conceptual challenges first:
What you have looked into, and already understand quite well, is the way in which FieldTrip does a lot of bookkeeping in order to try and ensure that all data objects that implicitly have a channel dimension are well matched w.r.t. the assumed order of the channels.
Specifically, in the context of forward and inverse modelling, different types of data ‘come together’. There’s leadfieds (either or not precomputed), spatial filters, and data (covariances and time series data). In order for the mathematical operations to make sense, the identity and the order of the channels in the data objects need to match. This assumed equivalence of channel order is not always guaranteed, particularly if you have a quirky naming (and ordering) scheme as the la Timone MEG system, but perhaps the other JM (=JeanMichel) would disagree with me calling your system quirky. Typically, the channel order in raw data from 4D-neuroimaging/Magnes/BTI systems is NOT in alphanumeric order, neither in the time series, nor in the description of the magnetometer array. Sooo, the long story short is that books need to be kept regarding the channel ordering. Yet, another thing to keep in mind in addition, is that any (linear) operations applied to the time series data (e.g. removal of spatial topographies by using ICA followed by backprojection, or applying denoising using the reference channels) need to be represented in the leadfields as well. In other words, if you precomputed the leadfields on an unbalanced grad-structure (fieldtrip speak for a sensor representation that has not been processed with the same processing pipeline that was used to preprocess the actual data) then this is not good.
Long story short, if you are to precompute your leadfields in a separate step (which I’d recommend), I would do it as follows:
cfg = [];
cfg.headmodel = <headmodel>
cfg.sourcemodel = <sourcemodel>
cfg. < some other options >
cfg.channel = data.label (assuming that the data does not contain non-MEG channels)
leadfield = ft_prepare_leadfield(cfg, data)
where data is the data structure that you later want to subject to ft_sourceanalysis. This would result in leadfields that not only have a more or less consistent order in the channels (i.e. consistent with the data), but also the leadfields are properly ‘balanced’. Note, however that a lot of the stuff you described in your e-mail related to the code that you identified in FT tries to ensure a correct channel order between different data objects.
Now, concretely, if you want to use a beamformer for source reconstruction, you should specify cfg.method = ‘lcmv’, rather than ‘mne’. If you on the other hand want to compute an mne-solution, you should NOT input the full data covariance matrix, but rather use an estimate of the noise (e.g. using an empty room measurement). In the code that you pasted in the e-mail you seem to have mixed these two things up.
Then, finally getting to what you want to achieve, the question is how to efficiently obtain single trial source level time series data. Indeed, I’d recommend a two-step approach, where you first obtain a set of spatial filters, and then subject the channel data to the spatial filters. I would assume that if you use precomputed leadfields as per the above suggestion, the order of the channels will match the order of the channels in the data (again provided the data does not contain non-MEG) channels. Now the question is whether you will have enough computer RAM to begin with to fit all the dipoles’ single trial time courses.
Perhaps we could follow up on this in a next set of e-mails in this thread - also to give other people the opportunity to chime in -, but also because 1) the above might already give you ample food for thought for now, and 2) my risotto is almost cooked, so I have to summon the troups at home to the dinner table.
Keep up the good work and critical thinking,
JM
On 25 Oct 2021, at 19:02, CHATEAUX Manon via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>> wrote:
Dear FieldTrip community,
I want to compute coherence between brain sources and an EMG channel that I selected. Here is the pipeline I am using :
* Preprocess MEG and EMG data
* Process single subject MRI data to get (a) the head model, (b) the source model with leadfields
* Reconstruct MEG data using beamformer techniques in the time domain.
* Compute coherence between source data and EMG data
To be able to compute coherence, I need to reconstruct data for all trials and not as an average of trials. But I encountered a problem at this level.
Here is the code I used first :
% Compute noise-covariance matrix
cfg = [];
cfg.covariance = 'yes'; %default takes the whole time window to estimate the noise covariance matrix
cfg.keeptrials = 'yes';
tlcontr_MEG = ft_timelockanalysis(cfg, MEG_contr);
% Compute inverse modelling with mne : source data averaged over
% trails
cfg = [];
cfg.method = 'mne';
cfg.sourcemodel = leadfield;
cfg.headmodel = headmodel;
cfg.mne.lambda = 3; %%????
cfg.mne.keepfilter = 'yes';
cfg.rawtrial = 'no';
source_MEG_contr = ft_sourceanalysis(cfg,tlcontr_MEG);
% Get source data for all trials : project trials by hand
FILTER = source_MEG_contr.avg.filter(cellfun(@(a) ~isempty(a),source_MEG_contr.avg.filter));
for triali = 1 : length(MEG_contr.trial)
source_MEG_contr_trials.pow{triali} = cell2mat(cellfun(@(a) sum((a*MEG_contr.trial{triali}).^2,1),FILTER,'UniformOutput',false)');
disp(['Trial n° ' num2str(triali) ' out of ' num2str(length(MEG_contr.trial))])
end
source_MEG_contr_trials is to be the structure that contains source data for all trials. As you can see I computed it by hand, following a thread that I found on the website.
Then I found it was possible to use ft_sourceanalysis again, to project all trials through the filter computed at the first use of the function. So I replaced the last paragraph of code with the following :
% project all trials through common spatial filter to get source data
% for separate trials
cfg=[];
cfg.method = 'mne';
cfg.sourcemodel = leadfield; % previously computed grid
cfg.headmodel = headmodel; % previously computed volume conduction model
cfg.rawtrial = 'yes'; % project each single trial through the filter to reconstruct single trial data
cfg.sourcemodel.filter = source_MEG_contr.avg.filter; % use the common filter computed in the previous step!
source_MEG_contr_trials = ft_sourceanalysis(cfg, tlcontr_MEG); % contains the source estimates for all trials/both conditions
I preferred this “automatic” way rather than the manual one because it allows me to keep the FieldTrip shape of structures more easily and display results using FieldTrip functions.
However, I realized that it returned completely different results. After trying to understand the steps performed in the second use of ft_sourceanalysis, I realized that filters used to compute trial source data from trial channel data (l. 295 in ft_inverse_mne of fieldtrip-20210921 : estimate.mom{i} = sourcemodel.filter{i}*dat) were different from filters in cfg.sourcemodel.filter.
In fact, channel labels are taken from the leadfield structure. Lines (channel dimension) in filter matrices (in cfg.sourcemodel.filter) are reorganized so that it matches channel order from data.label (= tlcontr_MEG.label).
See l.404-410 in ft_sourceanalysis of fieldtrip-20210921 :
[dum, chansel] = match_str(data.label, sourcemodel.label);
sourcemodel.label = sourcemodel.label(chansel);
for i=1:numel(sourcemodel.filter)
if ~isempty(sourcemodel.filter{i})
sourcemodel.filter{i} = sourcemodel.filter{i}(:, chansel);
end
end
Thus, if I understand it correctly, this step assumes that lines in filter matrices (which come from source_MEG_contr structure) are organized in the same way as channels in leadfields. However, this is not actually the case here because
1. source_MEG_contr.avg.label and data.label (= tlcontr_MEG.label) are equal.
2. I don’t understand why but leadfields.label is different from the original alphanumerical channel organization – A1 A2 A3 etc…- found in raw data (as well as in source_MEG_contr, and in tlcontr_MEG).
So I think that in my case, lines in filter matrices should not be reorganized (when they are not, I get the same results as by hand) but as leadfield structure is also used when I compute source_MEG_contr from tlcontr_MEG, I am not sure I can trust source_MEG_contr either.
I would truly appreciate if someone could proofread and pinpoint some problems in this reasoning to help me solve this annoying problem. In the meantime I will try to recompute the leadfield structure so that channel order is the same everywhere.
Thank you in advance for your help and the time you will spend reading through this message!
Manon Châteaux
phD Student at Aix-Marseille University
_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://urldefense.com/v3/__https://doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!sfGi3q4aOeMCoe7W0_ghPeD-viaovjeAzPnJRvx05cVOf7UZKHz_GIupVsUV7rEW0MzP5TGOqE2mnKs$
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20211030/d346eaa6/attachment.htm>
More information about the fieldtrip
mailing list