ENH: extract ROI time series

[JD-Zhu]

1. MNI305to152.py: utility to help read out anatomical labels from source localisation results.
2. traditional_source_analysis_and_ROIs.py: run source analysis and save the STCs, then extract source timecourses from predefined ROIs.

[larsoner]

a2009s is fine but nowadays I recommend HCPMMP1 like

mne.datasets.fetch_hcpmmp1_parcellation(subjects_dir=subjects_dir)

then here you can use the fine-grained (~360 labels) "HCPMMP1" or coarse (~44 labels) "HCPMMP1_combined" atlas

[JD-Zhu]

Thanks! I'll take a look into this :)

[JD-Zhu]

@larsoner Could HCPMMP1 be used with volumetric source results? The 'extract_label_time_course' method requires a string path to an atlas (e.g. ‘aparc.a2009s+aseg.mgz’) when working with volumetric source - the atlases are located in 'freesurfer/subjects/fsaverage/mri' and there doesn't seem to be a HCPMMP option available. Is this only used for surface source space? Thanks!

[larsoner]

@JD-Zhu no not directly, but some years ago I did write this

https://gist.github.com/larsoner/8e664205cd8285ca7c46211403ad12ce

the code comments at the top allow you to create a volumetric atlas that uses HCPMMP1 for the cortical regions and aseg for the subcortical regions, then that script creates an ID lookup that can be used with extract_label_time_course in volumetric mode. Want to take a look and give it a shot?

[JD-Zhu]

@larsoner Thanks, I'll have a go at it!
We were also hoping to run the ROI analysis with MEG-only for comparison - would the covariance matrix need to be recreated in this case (or could I just regenerate the fwd model and inverse solution based on MEG channels only)? If so, how do you make the cov matrix from the resting state data (I think this is what was done in the pipeline)?

[JD-Zhu]

Also quick question about the volumetric atlas for HCPMMP1 - there are some labels that sound very similar to each other, e.g.
'ctx-lh-precuneus' and 'ctx-lh-G_precuneus',
'ctx-lh-temporalpole' and 'ctx-lh-Pole_temporal'.
What is the difference between them and which ones should we use? Thanks so much!

[larsoner]

What is the difference between them and which ones should we use? Thanks so much!

For stuff like this I first try to check by loading the labels and doing brain.add_label, e.g.:

import mne
brain = mne.viz.Brain("fsaverage", hemi="lh")
labels = mne.read_labels_from_annot("fsaverage", "HCPMMP1")
labels = [label for label in labels if label.hemi == "lh" and "precuneus" in label.name]
for label in labels:
    brain.add_label(label, borders=True)

But this plots no labels! And checking we see len(labels) == 0, i.e, HCPMMP1 does not have labels by this name.

I suspect that you are seeing these names when processing FreeSurferColorLUT (right?) -- this is because that file has labels for a bunch of possible different existing FreeSurfer parcellations, some of which overlap. So for example looking at the two possible builting cortical parcellations you could use we see:

>>> mne.read_labels_from_annot("fsaverage", "aparc", regexp="precuneus", hemi="lh")
[<Label | fsaverage, 'precuneus-lh', lh : 7308 vertices>]
>>> mne.read_labels_from_annot("fsaverage", "aparc.a2005s", regexp="precuneus", hemi="lh")
[<Label | fsaverage, 'G_precuneus-lh', lh : 3558 vertices>]

And if you want to see where they are and/or how they overlap you can adapt the code above:

import mne
brain = mne.viz.Brain("fsaverage", hemi="lh", surf="inflated", views="medial")
label = mne.read_labels_from_annot("fsaverage", "aparc", regexp="precuneus", hemi="lh")
assert len(label) == 1, len(label)
brain.add_label(label[0], alpha=0.8)
label_a2005 = mne.read_labels_from_annot("fsaverage", "aparc.a2005s", regexp="precuneus", hemi="lh")
assert len(label_a2005) == 1, len(label_a2005)
brain.add_label(label_a2005[0], borders=True)

you can see that they are similar but the a2005 one sticks more to the gyral parts of the label (as implied by its G_ prefix).

[JD-Zhu]

Thanks! I was able to replicate your plot and could see how those two labels differ slightly. However, I'm still a bit confused as to how these labels from "aparc" and "aparc.a2005s" match with the labels in the .mgz file generated from your script linked above.

When I read in the .mgz file like this:
fname_aseg = subjects_dir / subject / 'mri' / 'HCPMMP1+aseg.mgz'
label_names = mne.get_volume_labels_from_aseg(fname_aseg)

I get a list of label names (strings, not actual Label objects). Is there a way to plot certain labels from this list to see where they are (similar to the plot you did above)? Sorry if I sound totally confused - still trying to grasp how these atlases and ROIs work!

[larsoner]

However, I'm still a bit confused as to how these labels from "aparc" and "aparc.a2005s"

So aparc, aparc.a2005s, HCPMMP1, HCPMMP1_combined, etc. are all cortical parcellations. These can be combined with aseg.mgz to create a volumetric atlas for the brain. freesurfer does this automatically for its built-in parcellations, e.g. aparc+aseg.mgz and aparc.a2005s+aseg.mgz are two volumetric atlases that are created for every recon-all'ed subject (including fsaverage) that correspend to two different cortical parcellations, each with aseg tacked on to handle subcortical regions.

So when you look at / think about HCPMMP1+aseg.mgz, you can view this in whatever MRI viewer you like using to look at volumetric atlases. One way is to use freeview the way you would any volumetric atlas. For aparc+aseg you can do:

freeview -v T1.mgz -v aparc+aseg.mgz:colormap=lut:opacity=1

For HCPMMP you'll probably need to pass a filename to the :lut option:

       ':lut=name' Set lookup table to the given name. Name can be the name of
       a stock color table or the filename of a color table file.

If you do this, you should see that the labels show up in the same places in the volumetric space that they do when you plot them on the inflated brain using FreeView or mne.viz.Brain(...) + brain.add_parcellation("HCPMMP1") for example.

[JD-Zhu]

Thanks for the explanation! I can now load the HCPMMP1+aseg.mgz into freeview with the LUT file from your script to see the locations.

Also, when I was reading in the label names previously, I was using the incorrect LUT file (the freesurfer default) so the label names were all incorrect which caused part of my confusion. I have now changed it to the following, which should load the labels correctly (I think):
fname_aseg = subjects_dir / subject / 'mri' / 'HCPMMP1+aseg.mgz'
fname_lut = subjects_dir / subject / 'mri' / 'HCPMMP1ColorLUT.txt'
lut = mne.read_freesurfer_lut(fname_lut)
label_names = mne.get_volume_labels_from_aseg(fname_aseg, atlas_ids=lut[0])

[JD-Zhu]

Hi @larsoner ,

I think I'm still having some trouble with the volumetric version of the HCPMMP1 atlas. The label names load correctly when using the proper LUT file, but I still run into a problem when trying to use these labels with mne.extract_label_time_course(). Apparently the label name cannot be found in the atlas even though it shows up in the list:

Does the custom lut file need to be supplied to mne.extract_label_time_course()? I checked the documentation but couldn't find a way to do so.

Thanks in advance for any advice.

[larsoner]

Does the custom lut file need to be supplied to mne.extract_label_time_course()? I checked the documentation but couldn't find a way to do so.

Not exactly, this one is a bit tricky -- unpacking the extract_label_time_course docs:

It looks like you have passed a two-element tuple, the secend element being a list of str. Following the hint to look at setup_volume_source_space, you'll see this note about labels:

You are in the "use of other atlases" part, so instead of passing a list of str you need to pass a dict whose keys are the volume label names and the values are integer IDs used in the atlas. You should be able to get those integer IDs from your custom LUT. Does that make sense?

[JD-Zhu]

You are in the "use of other atlases" part, so instead of passing a list of str you need to pass a dict whose keys are the volume label names and the values are integer IDs used in the atlas. You should be able to get those integer IDs from your custom LUT. Does that make sense?

All works now - thank you!