1.1 Reading in registered images

Here we will be reading in those previous registered, intensity-normalized images.

library(neurobase)
mods = c("T1", "T2", "FLAIR")
norm_reg_files = file.path("..", 
                     "preprocess_mri_within", 
                     paste0("113-01-", mods, "_norm_eve.nii.gz")
                     )
names(norm_reg_files) = mods
norm_reg_imgs = lapply(norm_reg_files, readnii)

1.2 Reading in Eve brain

Here we will read in the brain-extracted Eve T1 image, the brain mask, and then mask the normalized images with this mask.

library(EveTemplate)
eve_brain_fname = getEvePath("Brain")
eve_brain = readnii(eve_brain_fname)
eve_brain_mask = readEve(what = "Brain_Mask")
norm_reg_imgs = lapply(norm_reg_imgs, mask_img, mask = eve_brain_mask)

We will plot the registered subject images against this to ensure they are in fact in the same space. (Remember to always look at your data!)

lapply(norm_reg_imgs, double_ortho, x = eve_brain)

$T1
NULL

$T2
NULL

$FLAIR
NULL

We see good congruence from the template and the corresponding images from this patient.

1.3 Getting the Eve atlas and labels

Here we will read in one of the whole brain atlases of Eve, specifically the type 2 (“II”). Please see Oishi et al. (2009) and Oishi, Faria, and Mori (2010) for further discussion of the atlases. We will read in the atlas and show that it has a series of integer labels. These labels correspond to the integer_label column in the data.frame produced by getEveMapLabels.

eve_labels = readEveMap(type = "II")
unique_labs = eve_labels %>% 
  c %>% 
  unique %>% 
  sort 
head(unique_labs)

[1] 0 1 2 3 4 5

lab_df = getEveMapLabels(type = "II")
head(unique(lab_df$integer_label))

[1] 0 1 2 3 4 5

all( unique_labs %in% lab_df$integer_label)

[1] TRUE

1.4 Plotting the labels

Let’s plot the labels to see how they look.

ortho2(eve_labels)

Although that shows us the breakdown of labels, the defaults in ortho2 are for grayscale images. Luckily, getEveMapLabels provides an rgb representation of a look up table (LUT) for colors for each structure. Let’s use that to plot the image:

cols = rgb(lab_df$color_r/255, lab_df$color_g, 
           lab_df$color_b, maxColorValue = 255)
breaks = unique_labs
ortho2(eve_labels, col = cols, breaks = c(-1, breaks))

Again, not that great, but a bit better. We can try our own palette as below:

library(RColorBrewer)
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral')))
cols <- rf(length(unique_labs))
ortho2(eve_labels, col = cols, breaks = c(-1, breaks))

That may give us a good assortment of colors, where right is more red and blue is more left. We can randomize these colors, which may show better discrepancy by putting those with a blue hue close to those with red:

set.seed(20161008)
ortho2(eve_labels, col = sample(cols), breaks = c(-1, breaks))

We can overlay these colors on top of the subject’s T1 image:

set.seed(20161008)
ortho2(norm_reg_imgs$T1, eve_labels, 
       col.y = alpha(sample(cols), 0.5), ybreaks = c(-1, breaks))

1.5 Getting structure-specific metrics

Now that we have our normalized images in the template/atlas space, let’s get information about each imaging sequence. Here we will make a data.frame with all the voxels from each modality and the atlas. There are more elegant ways to do this, but we want to be explicit below:

df = data.frame(T1 = norm_reg_imgs$T1[ eve_brain_mask == 1],
                T2 = norm_reg_imgs$T2[ eve_brain_mask == 1],
                FLAIR = norm_reg_imgs$FLAIR[ eve_brain_mask == 1],
                integer_label = eve_labels[ eve_brain_mask == 1]
                )

Now that we have a standard data.frame, we can use any data manipulation procedure just as we would in any other way. Let’s reshape the data:

library(reshape2)
library(dplyr)
long = reshape2::melt(df, 
                      id.vars = "integer_label")
head(long)

  integer_label variable        value
1           130       T1  0.000000000
2           130       T1  0.000000000
3            65       T1 -0.006485715
4           130       T1  0.000000000
5           130       T1 -0.033087343
6           130       T1 -0.046731040

Now we can merge in the labels from the Eve template so that we can actually see what structures these voxels represent. We can drop the extraneous color columns as well.

long = left_join(long, lab_df, by = "integer_label")
long = long %>% 
  select(-color_r, -color_g, -color_b)

Now we can calculate some statistics:

stats = long %>% 
  group_by(integer_label, text_label, right_left, structure, variable) %>% 
  summarise(mean = mean(value),
            median = median(value),
            sd = sd(value)) %>% 
  select(variable, text_label, mean, median, sd, everything()) %>% 
  ungroup
head(stats)

# A tibble: 6 x 8
  variable                    text_label         mean      median
    <fctr>                         <chr>        <dbl>       <dbl>
1       T1                    background -0.606879821 -0.69270784
2       T2                    background  1.005639427  0.21436921
3    FLAIR                    background -0.660865609 -0.32353108
4       T1 superior_parietal_lobule_left  0.105027506  0.16597870
5       T2 superior_parietal_lobule_left  0.003699062 -0.38462870
6    FLAIR superior_parietal_lobule_left  0.037216124  0.03221552
# ... with 4 more variables: sd <dbl>, integer_label <dbl>,
#   right_left <chr>, structure <chr>

library(stringr)
stats %>% filter(str_detect(structure, "thalamus")) %>% 
  select(text_label, variable, mean) %>% 
  arrange(variable, text_label)

# A tibble: 6 × 3
      text_label variable        mean
           <chr>   <fctr>       <dbl>
1  thalamus_left       T1  0.80004855
2 thalamus_right       T1  0.78460175
3  thalamus_left       T2 -0.41903584
4 thalamus_right       T2 -0.39373081
5  thalamus_left    FLAIR  0.02892382
6 thalamus_right    FLAIR  0.05632016

transparent_legend =  theme(
  legend.background = element_rect(
    fill = "transparent"),
  legend.key = element_rect(fill =
                              "transparent", 
                            color = "transparent")
)

 long %>% filter(str_detect(structure, "thalamus")) %>% 
  ggplot(aes(x = value, colour = factor(right_left))) + 
  geom_line(stat = "density") + facet_wrap(~variable, ncol = 1) +
    theme(legend.position = c(0.2, 0.5),
        legend.direction = "horizontal",
        text = element_text(size = 24)) +
  guides(colour = guide_legend(title = NULL)) +
  transparent_legend

 long %>% filter(str_detect(structure, "thalamus")) %>% 
  ggplot(aes(x = variable, y = value, colour = factor(right_left))) + 
  geom_boxplot() +
    theme(legend.position = c(0.2, 0.85),
        legend.direction = "horizontal",
        text = element_text(size = 24)) +
  guides(colour = guide_legend(title = NULL)) +
  transparent_legend

2.1 Registration of Eve to the subject T1

Here we will use the N4 corrected, processed images from Processing Within-Visit MRI and register the Eve T1 brain to the processed T1 image. Let us get the filenames that we want:

eve_stub = "Eve_to_Subject"
outfile = paste0(eve_stub, "_T1.nii.gz")
lab_outfile = paste0(eve_stub, "_Labels.nii.gz")
outfiles = c(outfile, lab_outfile)

n4_files = file.path("..", 
                     "preprocess_mri_within", 
                     paste0("113-01-", mods, "_proc_N4_SS.nii.gz")
                     )
names(n4_files) = mods

Now we can register Eve to the subject T1 (note this is the opposite from Processing Within-Visit MRI). We will use a NearestNeighbor interpolator here. For the Eve brain image, this may not be a good interpolation and usually we would use a windowed sinc or linear interpolator.

We are using NearestNeighbor here for the labels, so that after interpolation, the value of the voxel in subject space is from the nearest neighbor, which ensures that the label is a integer within the origninal set. For example, let’s say after registration, a voxel is near labels 1 and 5. If we used an averaging interpolator (windowed sinc or linear), it may give the value of 3, which is an entirely different label and may not represent anything remotely close to correct. Moreover, these are labels so the numbering is arbirtrary, categorical values. Let’s say we shuffled the labels on the Eve atlas (but kept track of the structures), then the atlas is exactly the same. But if we applied the transformation to this shuffled atlas, the voxel may be near labels 420 and 6, which average to 120 or something. Thus, we want the interpolation for the labels to result in picking from a label set, which we generall would use NearestNeighbor or multiLabel, see antsApplyTransforms.

if ( !all( file.exists(outfiles) )) {
  reg = extrantsr::registration(
    filename = eve_brain_fname,
    template.file = n4_files["T1"],
    other.files = eve_labels,
    other.outfiles = lab_outfile,
    interpolator = "NearestNeighbor",
    typeofTransform = "SyN")
} 

# Running Registration of file to template

# Applying Registration output is

$warpedmovout
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 111x172x132 
  Voxel Spacing       : 1.20000004768372x1x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 


$warpedfixout
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 181x217x181 
  Voxel Spacing       : 1x1x1 
  Origin              : 0 0 0 
  Direction           : 1 0 0 0 -1 0 0 0 1 


$fwdtransforms
[1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpmTGkqd/file1323515f9ef4e1Warp.nii.gz"      
[2] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpmTGkqd/file1323515f9ef4e0GenericAffine.mat"

$invtransforms
[1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpmTGkqd/file1323515f9ef4e0GenericAffine.mat" 
[2] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpmTGkqd/file1323515f9ef4e1InverseWarp.nii.gz"

# Applying Transformations to file

# Applying Transforms to other.files

# Writing out file

# Writing out other.files

# Removing Warping images

# Reading data back into R

Note, the registration function is somewhat limited as only one interpolator can be specified. We could break the process into registration (with an averaging interpolator for the Eve brain) and then use extrantsr::ants_apply_transforms to apply this transformation to the Eve labels and use a nearest-neighbor interpolator. This is how it is set up in ANTsR and the restriction of the registration function comes at the expense of having one function that performs registration and interpolation in one. In this example, we don’t care specifically about the registered-to-subject Eve template, so we don’t break the process into two.

2.2 Overlaying Native Eve Labels on T1

Here we will read in the output from the registered Eve labels:

n4_proc_imgs = lapply(n4_files, readnii)
native_eve_labels = readnii(lab_outfile)

We will overlay the native Eve labels on top of the native-space N4 T1 images.

set.seed(20161008)
ortho2(n4_proc_imgs$T1, native_eve_labels, 
       col.y = alpha(sample(cols), 0.5), 
       ybreaks = c(-1, breaks))

2.3 Native data.frame

native_mask = n4_proc_imgs$T1 > 0
norm_imgs = lapply(n4_proc_imgs, zscore_img, 
                   mask = native_mask)
native_df = data.frame(
  T1 = norm_imgs$T1[ native_mask == 1 ], 
  T2 = norm_imgs$T2[ native_mask == 1 ],
  FLAIR = norm_imgs$FLAIR[ native_mask == 1 ],
  integer_label = native_eve_labels[ native_mask == 1 ]
)
nlong = reshape2::melt(native_df, 
                      id.vars = "integer_label")

Now we can merge in the labels from the Eve template so that we can actually see what structures these voxels represent. We can drop the extraneous color columns as well.

nlong = left_join(nlong, lab_df, 
                  by = "integer_label")
nlong = nlong %>% 
  select(-color_r, -color_g, -color_b)

g = nlong %>% filter(str_detect(structure, "thalamus")) %>% 
  ggplot(aes(x = value, colour = factor(right_left))) + 
  geom_line(stat = "density") + facet_wrap(~variable, ncol = 1) +
    theme(legend.position = c(0.2, 0.5),
        legend.direction = "horizontal",
        text = element_text(size = 24)) +
  guides(colour = guide_legend(title = NULL)) +
  transparent_legend
g + ggtitle("Native Space Thalamus")

(g + ggtitle("Eve Space Thalamus")) %+% 
  (nlong %>% filter(str_detect(structure, "thalamus")))

rm(list = "g")