Montage coordinate frames are saved wrongly

[moritz-gerster]

Dear MNE_BIDS Team,

Thank you for writing great software to easily handle data standardization! 😊 Everything works fine but our research group is struggling with one annoying bug that is hard to deal with:

Describe the bug

We save the montage for our raw data like this:

montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec)), 
                                                  coord_frame='mni_tal')
raw.set_montage(montage, on_missing='warn')
mne_bids.write_raw_bids(raw=raw, bids_path=bids_path, overwrite=True)

However, when we load the data, process it, and save it again, we get an error:

bids_path_load = BIDSPath(subject="003", session="EphysMedOn01", task="Rest", acquisition="StimOff", run="01",
                          processing=None, recording=None, space=None, split=None, root="../../rawdata",
                          suffix="ieeg", extension=None, datatype="ieeg", check=True)
raw = read_raw_bids(bids_path_load, verbose=False)
raw.resample(128)
bids_path_save = bids_path_load.copy().update(root="../../derivatives", processing="downsample")
write_raw_bids(raw, bids_path_save, allow_preload=True, format="BrainVision")

Expected results

Saving works:

Actual results

Saving does not work and yields an error message:

The reason this does not work is that the coordinate system is set to "head" even though we saved it as "mni_tal". If we set it to "mni_tal" again, it now works:

montage = raw.get_montage()
new_montage = mne.channels.make_dig_montage(ch_pos=montage.get_positions()["ch_pos"], coord_frame='mni_tal')
raw.set_montage(new_montage)
write_raw_bids(raw, bids_path_save, allow_preload=True, format="BrainVision")

How can we make sure we save the data correctly and avoid this reset of the montage?

Our coordinate system .json file looks like this:

sub-003_ses-EphysMedOn01_space-MNI152NLin2009bAsym_coordsystem.json:
{
"IntendedFor": "n/a",
"iEEGCoordinateProcessingDescription": "Co-registration, normalization and electrode localization done with Lead-DBS",
"iEEGCoordinateProcessingReference": "Horn, A., Li, N., Dembek, T. A., Kappel, A., Boulay, C., Ewert, S., et al. (2018). Lead-DBS v2: Towards a comprehensive pipeline for deep brain stimulation imaging. NeuroImage.",
"iEEGCoordinateSystem": "Other",
"iEEGCoordinateSystemDescription": "MNI152 2009b NLIN asymmetric T2 template",
"iEEGCoordinateUnits": "mm"
}

We would be really thankful if you could have a look at that! 👌🙂

Additional information

Platform: macOS-11.6.3-x86_64-i386-64bit
Python: 3.10.2 | packaged by conda-forge | (main, Feb 1 2022, 19:30:18) [Clang 11.1.0 ]
Executable: /opt/anaconda3/envs/bids_neu/bin/python
CPU: i386: 8 cores
Memory: Unavailable (requires "psutil" package)
mne: 0.24.1
numpy: 1.22.2 {blas=NO_ATLAS_INFO, lapack=lapack}
scipy: 1.8.0
matplotlib: 3.5.1 {backend=module://matplotlib_inline.backend_inline}

sklearn: Not found
numba: Not found
nibabel: 3.2.2
nilearn: Not found
dipy: Not found
cupy: Not found
pandas: 1.4.1
mayavi: Not found
pyvista: Not found
pyvistaqt: Not found
ipyvtklink: Not found
vtk: Not found
PyQt5: Not found
ipympl: Not found
mne_qt_browser: Not found
pooch: v1.6.0
mne_bids: 0.9

[welcome]

Hello! 👋 Thanks for opening your first issue here! ❤️ We will try to get back to you soon. 🚴🏽‍♂️

[hoechenberger]

Hello @moritz-gerster, thank you for reporting this!

@sappelhoff @adam2392 I suppose you two might have an idea?

[sappelhoff]

I know it's related to #723, but I don't know the fix right now. It'll probably take a relatively deep dive into the code. I could try to do that some time in April or May unless somebody beats me to it. (Sorry for the long delay)

PS: thanks for the detailed report Moritz!

[moritz-gerster]

Thanks! Yes, I am collaborating with Timon Merk who wrote #723 one year ago. It would be really great to overcome this issue at some point ☺️

[agramfort]

@moritz-gerster can you a standalone code snippet to replicate the problem? with what can be used as a unit test it will speed up the bug fix a lot. thx

[moritz-gerster]

Dear @agramfort, I expected you to raise this. We are working on an replication with some public mne-dataset!

[agramfort]

I am so predictable :) Message ID: ***@***.***>

…

[hoechenberger]

I am so predictable :)

We just have perfectly trained classifiers at our disposal!

[timonmerk]

Hello everyone, thanks for picking up on this issue! The PR I created reads with raw_read.get_montage().get_positions()["coord_frame"] to "other" instead of "MNI152NLin2009bAsym", so the assert fails.

[moritz-gerster]

Dear @agramfort, is this PR sufficient for pinpointing the issue? :)

[alexrockhill]

Hey @moritz-gerster, sorry in being late to the party on this issue, this is something that I changed recently so I think I'll have your answer.

But first, a few questions: which coordinate frame are your ieeg channels in? Do you want to save them in a template (mni_tal) or subject (mri i.e. in relation to the subject's anatomy)?

For some context, the groups that developed the BIDS standard were using "other" as their coordinate frame to mean that it was in subject MRI space (https://bids-specification.readthedocs.io/en/stable/99-appendices/08-coordinate-systems.html#ieeg-specific-coordinate-systems). This was very confusing to me and seemed to defeat the whole purpose of BIDS so we made some clarifications and put them in this tutorial (https://mne.tools/mne-bids/dev/auto_examples/convert_ieeg_to_bids.html). If you are making a minimally reproducible example at some point, it could be helpful to hack the tutorial (link just before) so that it matches your data and then I could replicate on my machine.

[alexrockhill]

Ah I read more carefully and you're trying to use MNI152NLin2009bAsym which is a template space. You should just be able to follow the directions here https://mne.tools/mne-bids/dev/auto_examples/convert_ieeg_to_bids.html#step-5-store-coordinates-in-a-template-space (by using apply_trans on your montage to make sure it's in the mni_tal coordinate frame).

I think the only limitation is that the dataset will not say that the coordinates are in terms of MNI152NLin2009bAsym but you could go in and change that by hand in the coordinate_system.json file. We should add an argument to do this programmatically though.