AAL2 and AAL3 atlas

[drmowinckels]

Hello Dr. Mowinckel,

Is the AAL2 atlas available in ggseg or are there plans to include this atlas? Thanks!

Best,
John

Originally posted by @neuropil in ggseg/ggsegExtra#12 (comment)

[drmowinckels]

https://www.gin.cnrs.fr/en/tools/aal/

[fcmeyer]

Hi Dr. Mowinckel (@Athanasiamo),

Firstly, let me thank you for this AMAZING software you are developing! As an avid userR and neuroimager, finding ggseg was like a dream come true 😍 And it was exciting to see how many atlases you have added in such a short span of time!!

So, I have been trying to contribute to your project by building the AAL3v1 atlas, which is what I have been using for my dissertation. Unfortunately, I have gotten stuck and stumped in this process. I am not a big freesurfer person, so perhaps I am doing something wrong here...

To give you some background, this atlas is one of those trickier cases (from what I gathered from perusing your code / github discussions) because (1) only provided in volumetric; (2) has both cortical and subcortical ROIs. I have done my best to approach this in numerous ways and wanted to share where I am so far to see if you may have any suggestions on moving forward.

Existing code & issues

Below, I share some code I created for this purpose. What works and doesn't:

• Generating a ggseg3d atlas using make_volumetric_2_3datlas()
• Plot 3d atlas using ggseg3d
• Plotting data on 3d atlas
• Separating subcortical from cortical, and hemispheres
• Attempting to generate a ggseg 2d atlas from 3d using make_ggseg3d_2_ggseg()
• Attempting to generate a ggseg 2d atlas using make_volumetric_ggseg()

Please read on below this code chunk, where I will elaborate on some of these issues.

library(ggsegExtra) # obvious reasons
library(RColorBrewer) # for inventing a color palette
library(magrittr) # for piping 
require(here) # for paths
require(stringr) # for various things
require(tidyr)
require(dplyr)
require(tibble)
require(glue)

# The unique name of the atlas annot, without hemisphere in filename
annot_name <- "aal3v1"

# Download AAL3v1 and extract --------------------------------------------------
url <- 'https://www.oxcns.org/AAL3v1_for_SPM12.zip'
tmp <- tempfile()
download.file(url, tmp)
unzip(tmp, exdir=here::here('data-raw'))

# Register atlas to fsaverage5 / MNI305 space ----------------------------------

# Filenames for Freesurfer-compatible files to be generated
label_file <- here::here('data-raw','AAL3v1_fs5.mgz')
lut_file <- here::here('data-raw', 'AAL3v1_lut.txt')

# The template is provided in various formats; according to the docs, the
# ROI_MNI_V7 and AAL3v1 files have the same data, jsut different headers and 
# compression; here I am using the uncompressed NIFTI file AAL3v1.nii
vol2vol_source_file <- here::here('data-raw','AAL3','AAL3v1.nii')

# Command to run in bash to create file in MGZ format / fsaverage5 space
# Using line breaks for readability :)
fs_cmd <- "mri_vol2vol
--reg $FREESURFER_HOME/average/mni152.register.dat
--nearest
--mov {vol2vol_source_file}
--targ $FREESURFER_HOME/subjects/fsaverage5/mri/brain.mgz
--o {label_file}"

# Clean up line breaks for use in system, fill in paths
fs_cmd <- fs_cmd %>% 
  stringr::str_replace_all('\n',' ') %>% 
  glue::glue()

# For my Mac only, because my paths are weird...
# options(freesurfer.path='/Applications/freesurfer/7.2.0/')
# withr::with_envvar(
#   new = c(FREESURFER_HOME=getOption('freesurfer.path'),
#           SUBJECTS_DIR=file.path(getOption('freesurfer.path'),'subjects')),
#   system(glue::glue('source $FREESURFER_HOME/SetUpFreeSurfer.sh; {fs_cmd}'))
# )

# For computers where freesurfer is correctly in your path and you don't have
# all this drama...
system(fs_cmd)

# Built LUT -------------------------------------------------------------------

# Read labels and build LUT
lut <- vroom::vroom(here::here('data-raw','AAL3v1.nii.txt'),
                    col_names=FALSE)
names(lut) = c('idx','label','old_aal3_label')

# My hackish solution to invent a color palette... perhaps there is something
# better... let me know if so!
coul <- brewer.pal(11,'Spectral')
coul <- colorRampPalette(coul)

# Do table
lut <- lut %>%
  dplyr::select(-old_aal3_label) %>%
  dplyr::mutate(
    col_hex = coul(n())
  ) %>%
  dplyr::mutate(
    col = purrr::map(col_hex, ~ col2rgb(.x, alpha=TRUE) %>% 
                       t() %>% 
                       data.frame() %>% 
                       tibble::tibble()
    )
  ) %>%
  tidyr::unnest(cols=c('col')) %>%
  dplyr::select(-col_hex) %>%
  dplyr::rename(R=red, G=green, B=blue, A=alpha)
  
# Save LUT
write_ctab(lut, lut_file)

# Compute 3d Atlas -------------------------------------------------------------

# Compute atlas object
aal3v1_3d <- make_volumetric_2_3datlas(
  template = label_file,
  color_lut = lut_file,
  subject='fsaverage5',
  output_dir=here::here('data-raw'),
  ncores=8,
  cleanup=FALSE,
  verbose=TRUE
)

# Post processing steps --------------------------------------------------------

# Make palette and do the other things in the provided script
brain_pals <- make_palette_ggseg(aal3v1_3d)
usethis::use_data(brain_pals, internal = TRUE, overwrite = TRUE)
devtools::load_all(".")

usethis::use_data(aal3v1_3d,
                  internal = FALSE,
                  overwrite = TRUE,
                  compress="xz")

# Testing plotting -------------------------------------------------------------

# This works
ggseg3d(atlas = aal3v1_3d)

# But plotting data don't work...

# Making up data based on the ggseg3d vignettes
mydata <- aal3v1_3d %>%
  unnest(ggseg_3d) %>% 
  ungroup() %>% 
  select(region) %>% 
  na.omit() %>% 
  mutate(p = runif(n(),4,8))

# Attempt to plot made up data on the brain
ggseg3d(data = mydata, 
        atlas = aal3v1_3d,
        color = "p", text = "p")

Issues

Plotting data using ggseg3d (my highest priority, hehe)

Attempting to plot data using the last bit of the code chunk above yields the following error:

Error: Atlas must have surfaceinflatedLCBC

I am not sure if this might have to do with some freesurfer-related things that need to be done, or what... any advice is much appreciated!

Separating cortical from subcortical / hemispheres

If you look at the generated atlas, out of the box it shows as subcort hemisphere:

> aal3v1_3d
# A tibble: 1 × 4
  atlas         surf  hemi              ggseg_3d
  <chr>         <chr> <chr>   <list<tibble[,5]>>
1 AAL3v1_fs5_3d LCBC  subcort          [166 × 5]

I am not fully sure how to separate between what should be in subcort vs. lh/rh (sorry, freesurfer just ain't my expertise, haha), but my first (hackish and perhaps lazy) pass at it is as follows:

aal3v1_n <- aal3v1_3d
aal3v1_n <- unnest(aal3v1_n, ggseg_3d)
aal3v1_n <- mutate(aal3v1_n,
                    region = gsub("_L$|_R$", "", region),
                    hemi = ifelse(str_detect(label, '_L'),'left',
                                  ifelse(str_detect(label,'_R'),'right','subcort')),
                    atlas = "aal3v1_3d"
)
aal3v1_3d_lateralized <- as_ggseg3d_atlas(aal3v1_n)
ggseg3d(atlas  = aal3v1_3d_lateralized)

If you have any suggestions on more programatic to identify specifically "which" regions should go as lh/rh vs subcort, that would be much appreciated!

Attempting to generate a ggseg 2d atlas from 3d using make_ggseg3d_2_ggseg()

When running the following:

aal3v1 <- make_ggseg3d_2_ggseg(aal3v1_3d_for2d, 
                                 output_dir = here::here("data-raw/"))

I get this error:

Error: Atlas must have surfaceinflatedLCBC

I'm not sure exactly what this might refer to, I imagine more freesurfer stuff needs to be done? Any advice would be much appreciated...

Attempting to generate a ggseg 2d atlas using make_volumetric_ggseg()

At first, I had trouble with the default slice values. I ended up putting some random small values and that fixed the problem, for the time being. I imagine I will need to look at the coordinates and find the set of slices "with brain" to define this. Will give more thought to this later on.

Where I got fully stumped is Step 5. Here is the code I used (note: you probably want to run the code provided above up to the make_volumetric_to_3datlas() call to have all the necessary files and vars relied on here)

aal3v1 <- make_volumetric_ggseg(label_file = label_file
                                , subject='fsaverage5'
                                , color_lut=lut_file
                                ,slices = data.frame(x = c(0, 2, 3), 
                                                     y = c(0, 2, 3), 
                                                     z = c(0, 2, 3),
                                                     view = 
                                                       c("axial", "sagittal",
                                                         "coronal"), 
                                                     stringsAsFactors = FALSE)
                                , output_dir = here::here("data-raw/"),
                                , steps=c(1:8))

The program runs through step 5, and shows the progress bar, but then this happens...

%%  5/8  Extracting contours from regions
  |================================================================| 100%, Elapsed 08:07
Error in `group_by()`:
! Must group by variables found in `.data`.
x Column `filenm` is not found.

I wonder if this might be a bug in the code, or perhaps because I am doing something wrong. In my digging through your codebase, I noticed that filenm is defined for this pipeline in step 7, and wasn't able to find a place where this would occur prior for this pipeline. But perhaps it has to do with something odd with my data, who knows. In any case, any suggestions would be immensely appreciated!

Once again, thank you for your amazing work on this package, and also for taking the time to look through this issue. I look forward to hearing from you on any suggestions on how to move forward, so that I can hopefully contribute and share this with other researchers who would like to plot their ROI data for this atlas! 😄

Cheers,
Francisco Meyer

[fcmeyer]

I'm just realizing, I forgot to clarify an important thing - all the code I posted assumes you are using RStudio with a new ggsegExtra project loaded, hence the data-raw/ subdirectory assumption and use of here::here() calls. So, you would need to create a new ggsegExtra project on RStudio and then run the code from there.

Thanks again! :)

[drmowinckels]

this is great! This is the first proper outside debugging og atlas creation we have had. I'd love to keep testing and figuring things out.

Could you by chanse provide me with the atlas file, so I can try reproducing your work? Its the only file I'm missing I think to try doing it my self too. This should expose if I have something on my system making it work, or if there is something in the code bugging. I have a small suspicion on where the error might be, but cant really tell untill I give it a try locally.

[drmowinckels]

I take it back, I just saw the aatlas. url. i'll be doing a test :D Will keep you posted!

[drmowinckels]

Hi!

I'm making some progress here, but need to get some freesurfer issues sorted out (I have not used it since updating to monterey, and I cant even open freeview currently). So I cant get the pipeline to create a ggseg-atlas finalized. But the 3d version I could do some progress on.

As you noted, you do actually manage to make a 3d atlas, but it's not behaving exactly as you wanted. There are good reasons for this, though you might not like them 😛
When using the volumetric functions, which are necessary for atlases that include subcortical structures, there is no longer a real split between hemispheres. This is because there are subcortical structures spanning both hemispheres, and we made this decision to deal with atlases including subcortcal data in this way.

This is meaningful when you have pure subcortical atlases like the aseg atlas, but when you have a full brain atlas like this, it starts being annoying. I'll need to have a think about how I want to work with this. This is in fact the first full brain atlas ever attempted for ggseg to our knowledge, so its time to start thinking about it 👍🏼

When it comes to plotting your own data, there are a couple of misspells in your ggseg3d arguments stopping it from working. I maybe missed sections of the documentation with this when I changed the arguments, so I'm sorry if I made things confusing.

It should be .data and colour for the arguments. and then it works for me.

ggseg3d(.data = mydata, 
        atlas = aal3v1_3d,
        colour = "p", text = "p")

make_ggseg3d_2_ggseg will only work with a fully inflated surface and only on cortical atlases, which is why we made its companion make_volumetric_ggseg for atlases with subcortical data. But I cant figure this one out before I get freesurfer sorted, hopefully very soon.

[nicjohnware]

Good evening, team. Have you ever managed to get this sorted? I have just requested the AAL3 and stumbled upon this.

This would save us a lot of time! :)