cookbook_recipes.Rmd

---
title: "A recipe book of brain plots"
author: |
  | Yohan Yee
  | version 0.1
date: "March 4, 2021"
output: 
  rmdformats::readthedown:
    self_contained: true
    highlight: tango
    lightbox: true
    df_print: paged
urlcolor: blue
---

```{r setup, include=FALSE}
cwd <- "~/dev/miceviz/"

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(root.dir=cwd)
setwd(cwd)

library(tidyverse)
library(glue)
library(RMINC)
library(scales)
library(ggnewscale)
library(patchwork)
```

\newpage
# Introduction

This cookbook contains recipes to create various plots using `ggplot2`. Each recipe contains code that can be run independent of other recipes / cookbooks, and should re-create the featured plot given appropriate paths. 

## Get the ingredients

### Install the required packages

You will need the following packages:

* tidyverse
* glue
* RMINC
* scales
* patchwork
* ggnewscale

Note that `RMINC` requires the minc-toolkit, and Windows isn't supported. Consider using the [MINC virtual machine](https://github.com/CoBrALab/MINC-VM) (which comes with RMINC preintalled) if you are using Windows.

You can install packages using either `install.packages()`, or if the source code is on Github, via `install_github` from the `devtools` package:

```{r eval=FALSE}
# tidyverse, glue, and scales are available on CRAN
install.packages("tidyverse")
install.packages("glue")
install.packages("scales")
 
# Other packages are available on Github
# install.packages("devtools") # Uncomment if you don't have the devtools package and need to run the commands below
devtools::install_github("thomasp85/patchwork")
devtools::install_github("eliocamp/ggnewscale")
devtools::install_github("Mouse-Imaging-Centre/RMINC", upgrade_dependencies = FALSE)
```

### Download the cookbook ingredients

Data will be provided.

## ggplot2

`ggplot2` is a library for creating 2D plots. Building on the principle that visualizations can be expressed through a "grammar of graphics" (hence the name), `ggplot2` provides a language to build beautiful plots layer by layer, and after some use, you will find this language to be quite intuitive. 

A handy cheatsheet can be found [on GitHub](https://raw.githubusercontent.com/rstudio/cheatsheets/master/data-visualization-2.1.pdf), or if you are using RStudio, by clicking on *Help* > *Cheatsheets* > *Data Visualization with ggplot2*.

The subsections below are a ten minute introduction to plotting data with `ggplot2`.

### Basic plotting, grouping, and positioning

Load the data.

```{r}
# Output of anatGetAll() bound with subject metadata
# From Lindenmaier et al. (2021)
# Thanks Zsu for providing this
# Ignore duplicated column warning for now
df <- read_csv("resources/Lindenmaier_2021_NeuroImage_AndR_absvols.csv")

# Example of getting whole brain volumes using dplyr
whole_brain_volume <- df %>% 
  select(-Project, -ID, -StrainW, -Genotype, -Sex, -Unique) %>%
  rowSums()

df <- df %>%
  mutate(whole_brain_volume=whole_brain_volume) %>%
  select(Project:Unique, whole_brain_volume, everything())

print(df)
```

### Quick plotting

`ggplot` takes in a data frame and plots data points (rows) as geometric objects, using the `aes()` function to map columns to each part of the plot. Let's look at the distribution of whole brain volumes in this dataset.

```{r}

# Start by calling ggplot and mapping whole brain volume to the y-axis
# This will plot nothing for now
ggplot(data = df, mapping = aes(y=whole_brain_volume))

# We need to tell ggplot to add points on top
# Add geom_point(). geom_point() plots points at an x,y location and thus requires both
ggplot(data = df, mapping = aes(y=whole_brain_volume)) + geom_point(x=0.5)

# Location can be factor, rather than continuous
# Here, let's use genotype
# Since we're getting this information from the data frame (i.e. the genotype depends on the point / mouse, and is different between rows in the data frame),
#   we need to wrap it in aes()
ggplot(data = df, mapping = aes(y=whole_brain_volume)) + geom_point(mapping = aes(x=Genotype))

# geom_point() / additional layers inherit data from the ggplot() object if nothing specified, this can be overwridden
ggplot(data = df, mapping = aes(x=Genotype, y=whole_brain_volume)) + geom_point()
ggplot(data = df, mapping = aes(x=Genotype, y=whole_brain_volume)) + geom_point(mapping = aes(x=Sex))

# Similarly, data can also be manually specified
# Sample only 5 points and plot
# Sampling can be helpful when prototyping plots with many points (e.g. millions of voxels in the brain)
ggplot(data = df, mapping = aes(x=Genotype, y=whole_brain_volume)) + geom_point(data = df %>% sample_n(5))

# Data can also be piped into ggplot, and code organized between lines
# so that the code is readable
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) + 
  geom_point()

# And the plot can be made readable too, by jittering points
# In this case, we only want to jitter in the x (width)-direction to avoid overlaps, but not in the y (height)-direction, so that
#   the data are not changed
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) + 
  geom_point(position=position_jitter(width = 0.2, height=0))

# Jittering is random and changes if called again
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) + 
  geom_point(position=position_jitter(width = 0.2, height=0))

# We can also further map variables to colour
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume, color=Sex)) + 
  geom_point(position=position_jitter(width = 0.2, height=0))

# And dodge the points
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume, color=Sex)) + 
  geom_point(position=position_jitterdodge(jitter.width = 0.2, jitter.height=0, dodge.width = 0.5))

```

### Layers

There exist many other ways to plot these distributions by changing very few things:

```{r}

# Boxplot
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) +
  geom_boxplot()

# Dotplot
# Bins of 5 units in the y-axis
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) +
  geom_dotplot(binaxis = "y", stackdir = "center", binwidth = 5)

# Violin plot
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) +
  geom_violin()

# Violin plot but with points, using a function from the ggforce package
df %>%
  ggplot(mapping = aes(x=Genotype, y=whole_brain_volume)) +
  ggforce::geom_sina()

# Bar plot
df %>%
  group_by(Genotype) %>%
  summarize(mean_wbv=mean(whole_brain_volume)) %>%
  ungroup %>%
  ggplot(mapping = aes(x=Genotype, y=mean_wbv)) +
  geom_bar(stat = 'identity')

# Histogram
df %>%
  ggplot(mapping = aes(x=whole_brain_volume, fill=Genotype)) +
  geom_histogram()
  
# Density
df %>%
  ggplot(mapping = aes(x=whole_brain_volume, fill=Genotype)) +
  geom_density(alpha=0.5)
  
```

### Quick plotting

Let's look at how the volume of the midbrain scales with whole brain volume. Normally you would plot

```{r}
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point()
```

A quick interface to plotting points is the `qplot()` function:

```{r}
qplot(x=whole_brain_volume, y=midbrain, data = df)
```

### Facets

Data can be split across multiple panels. 

```{r}
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point() +
  facet_wrap(~ Sex)
```

Suppose we want to show multiple panels but for different structures (columns) in the dataset. How do we facet by columns? Answer: `gather()` the data into "long" form. To keep things simple, let's do this for four structures: left and right thalamus, midbrain, and pons

```{r}
# Only select columns of interest
df_selected <- df %>%
  select(Project:whole_brain_volume, ends_with("_thalamus"), midbrain, pons)

# Gather into long form
df_selected_long <- df_selected %>%
  gather(key = "Structure", value = "Volume", ends_with("_thalamus"), midbrain, pons)

# Facet by Structure now
# Set scales="free" otherwise the same x and y axes will be used for all plots, which isn't wanted since the structures differ in volume
df_selected_long %>%
  ggplot(aes(y=Volume, x=whole_brain_volume)) +
  geom_point() +
  facet_wrap(~ Structure, scales = "free")

# Can facet by multiple variables
df_selected_long %>%
  ggplot(aes(y=Volume, x=whole_brain_volume)) +
  geom_point() +
  facet_wrap(~ Structure + Sex, scales = "free")

# And organize into a grid
df_selected_long %>%
  ggplot(aes(y=Volume, x=whole_brain_volume)) +
  geom_point() +
  facet_wrap(~ Structure + Sex, scales = "free", nrow=4)

# Or facet by two variables, in a grid:
df_selected_long %>%
  ggplot(aes(y=Volume, x=whole_brain_volume)) +
  geom_point() +
  facet_grid(Sex ~ Structure)

```


### Stacking layers

We can stack multiple layers:

```{r}
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point() +
  geom_smooth(method='lm', formula = y ~ x) +
  facet_wrap(~ Sex)
```

And layer order matters.

```{r}
# Bring points to the front of the line
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_smooth(method='lm', formula = y ~ x) +
  geom_point() +
  facet_wrap(~ Sex)
```

### Scales and coordinates

Scales can be modified

```{r}
# Reverse x axis scale
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_smooth(method='lm', formula = y ~ x) +
  geom_point() +
  facet_wrap(~ Sex) +
  scale_x_reverse()

# Manual colour scale while keeping reversed x scale
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_smooth(method='lm', formula = y ~ x) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_reverse() +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'))
```

### Labels, titles, and legends

Many more things in the plot can be modified:

```{r}
# Axis breaks
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'))
  
# Axis titles
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20), name="Whole brain volume") +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5), name="Midbrain") +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'))

# Axis limits
# Note the warning when points are outside these limits
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20), name="Whole brain volume") +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5), name="Midbrain volume", limits=c(12,15)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'))
  
# Can also specify labels and limits outside scales
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40')) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  ylim(12, 15) 

# Plot title
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40')) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag")

# Legend title
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title") +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag")

# Remove legend
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title", guide=F) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag")
  
```


### Themes

Various themes can be used:

```{r}
# Dark
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title", guide=F) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_dark()

# Black and white
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title", guide=F) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_bw()

# Classic
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title", guide=F) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_classic()

# Linedraw
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title", guide=F) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_linedraw()

# Void
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title", guide=F) +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_void()
```

Theme elements can be manually specified:

```{r}
# Note: must specify manual elements after adding a theme, otherwise that theme will override manual specifications
df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title") +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_bw() +
  theme(axis.title = element_text(size=14),
        axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "bottom")
```

### Saving plots

Plots can be saved to disk. For publication plots, use 90 dpi or more. 

```{r}
# First set as variable
plt <- df %>%
  ggplot(aes(y=midbrain, x=whole_brain_volume)) +
  geom_point(aes(color=Genotype)) +
  facet_wrap(~ Sex) +
  scale_x_continuous(breaks=seq(from=400, to=520, by=20)) +
  scale_y_continuous(breaks=seq(from=11, to=17, by=0.5)) +
  scale_color_manual(values=c('#FFAAAA', 'red', 'grey40'), name="Manual\nlegend title") +
  xlab("Whole brain volume") +
  ylab("Midbrain volume") + 
  labs(title = "Title", subtitle = "Subtitle", caption = "Caption", tag = "Tag") +
  theme_bw() +
  theme(axis.title = element_text(size=14),
        axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "bottom")

# Save via ggsave
ggsave(plot = plt, filename = "example_plot_ggsave.png", device = png(), dpi = 90, width=6, height = 4, units = "in")

# Save by printing to device
png(filename = "example_plot_pngdevice.png", width = 2000, height = 1000, res = 150)
print(plt)
dev.off()

```

### Ploting a brain slice

<div class="alert alert-info">
  <strong>Note:</strong> the `MRIcrotome` package provides an easier interface for plotting brain slices. Plots can be customized with `ggplot` however.
</div>

```{r}

# Get files
DSURQE_average <- glue("resources/DSURQE_40micron_average.mnc")
DSURQE_labels <- glue("resources/DSURQE_40micron_labels.mnc")

# Read in data as an array
bgarray <- mincArray(mincGetVolume(DSURQE_average))[80:160,204,120:180] 
atlasarray <- mincArray(mincGetVolume(DSURQE_labels))[80:160,204,120:180]

# Set row and column names
bgarray <- bgarray %>% `dimnames<-`(list(as.character(80:160), as.character(120:180)))
atlasarray <- atlasarray %>% `dimnames<-`(list(as.character(80:160), as.character(120:180)))

# Convert data to long form for plotting
bgdf <- bgarray %>%
  as_tibble %>%
  mutate(x=rownames(bgarray)) %>%
  gather(z, intensity, -x) %>%
  mutate(x=as.numeric(x),
         z=as.numeric(z))

atlasdf <- atlasarray %>%
  as_tibble %>%
  mutate(x=rownames(atlasarray)) %>%
  gather(z, label, -x) %>%
  mutate(x=as.numeric(x),
         z=as.numeric(z))

# Plot
bgdf %>%
  ggplot(aes(x=x, y=z)) +
  geom_raster(aes(fill=intensity), interpolate = T) +
  scale_fill_gradient(low="black", high="white", limits=c(800, 1800), oob=squish, guide=F) +
  new_scale_fill() +
  geom_raster(aes(fill=factor(label)), 
              alpha=0.2,
              interpolate = F, 
              data=atlasdf %>%
                filter(label %in% c(327, 329, 331))) +
  scale_fill_brewer(palette = "Spectral", name="Label", labels=c("MoDG", "GrDG", "PoDG")) +
  coord_fixed() + 
  xlab("x-coordinate (voxel index)") +
  ylab("z-coordinate (voxel index)") +
  labs(title="Anatomy",
       subtitle = "with a bunch of hippocampal labels") +
  theme_void()

```




\newpage
# Brain slices on a white background

## Goal

Plot a series of coronal slices through the brain similar to `MRIcrotome`'s `sliceSeries()` or `RMINC`'s `mincPlotSliceSeries()` functions, but on a white background.

## Ingredients

Libraries:

```{r}
library(tidyverse)
library(glue)
library(RMINC)
library(scales)
library(patchwork)
```

Input arguments:

* `background_file`, a path to a .mnc file, corresponding to the background anatomy
* `mask_file`, a path to a .mnc file, that defines a brain mask corresponding to the anatomy
* `slices`, a vector of indices, corresponding to slices (based on world coordinates) that you want to plot

```{r}
background_file <- "resources/brain_slices/average_template_50.mnc"
mask_file <- "resources/brain_slices/average_template_50_mask.mnc"
slices <- seq(from=-3, to=2, by=1)
```

## Step 1: Load the data

Load the data (.mnc files corresponding to brain anatomy and mask) that you want to plot using `mincGetVolume()`. 

```{r}
# Read in the MINC volume 
vol <- mincGetVolume(background_file)
mvol <- mincGetVolume(mask_file) 
```

## Step 2: Fix non-brain values so that they will appear white

First, determine an intensity value that will be considered white. Then, set values outside values outside the brain mask to this value.

```{r}
# Set the background to white
# Idea: the background (stuff outside the brain mask) can be set to the highest value within the brain (white on a greyscale colour bar with standard range)
# We'll improve the idea by setting it to a value bit less than that, to improve contrast
max_intensity <- max(vol[mvol > 0.5])
vol[mvol < 0.5] <- 0.75*max_intensity

# Convert to 3D array
arr <- mincArray(vol)
```

## Step 3: Convert world coordinates to voxel indices

For each slice, convert world coordinates to voxel indices for the slices to be plotted. The output is a data frame with three columns corresponding to each slice's world and voxel coordinates, and an index. The data frame is printed out at the end to see this. 

```{r}
# Convert to voxel coordinates by mapping over world coordinates
# This is specifically for coronal slices
# Note: if modifying this code for other dimensions (e.g. sagittal slices), confirm the ordering of these dimensions
#   - some tools require/output coordinates as x,y,z while others use z,y,x ordering
#   - for mincConvertWorldToVoxel, the input is x, y, z (relative to Display), and returns voxel coordinates as z, y, x (zero indexed)
slices_voxel <- slices %>%
  map_dfr(
    function(w) {
      world_coord_y <- w
      voxel_coord_y <- mincConvertWorldToVoxel(background_file, 0, w, 0)[2]
      out <- tibble(world=world_coord_y, voxel=voxel_coord_y)
      return(out)
    }
  ) %>%
  arrange(desc(world)) %>%
  mutate(index=1:nrow(.))

# Print 
print(slices_voxel)
```

## Step 4. Get starting and ending (world) coordinates

Get starting and ending (world) coordinates so we can appropriately label the axes of the plot. 

```{r}
starts <- mincConvertVoxelToWorld(background_file, 0, 0, 0)
ends <- mincConvertVoxelToWorld(background_file, dim(arr)[3]-1, dim(arr)[2]-1, dim(arr)[1]-1)
```

## Step 5: Extract slices and put into data frame

Extract slices for these voxel coordinates and wrangle the data to long format. Slicing the array provides a 2D array of intensity values for each slice. To plot this with `ggplot` however, we need to put the data into "long" form, so that we have two columns corresponding to the row and column index, and a third column corresponding to intensity.

```{r}
# Map over each coronal slice
plt_df <- seq_along(slices) %>%
  map_dfr(
    function(i) {
      
      # Get world and voxel coordinates corresponding to slice
      w <- slices_voxel$world[i]
      s <- slices_voxel$voxel[i]
      index <- slices_voxel$index[i]
      
      # Get 2D array corresponding to background intensities
      coronal_slice_array <- arr[,(s+1),]
      
      # Label row and column names for this array according to the world coordinates
      rownames(coronal_slice_array) <- as.character(seq(from=starts[1], to=ends[1], length.out = dim(arr)[1]))
      colnames(coronal_slice_array) <- as.character(seq(from=starts[3], to=ends[3], length.out = dim(arr)[3]))
      
      # Wrangle data to long form
      out <- coronal_slice_array %>%
        as_tibble %>%
        mutate(x=rownames(coronal_slice_array)) %>%
        gather(key = "y", value="intensity", -one_of("x")) %>%
        mutate(x=as.numeric(x), 
               y=as.numeric(y),
               coronal_slice_voxel=s,
               coronal_slice_world=w,
               coronal_slice_index=index)
      
      # Return 
      return(out)
    }
  )
```

## Step 6: Create a basic plot

Using the `geom_raster()` layer, we can plot the intensities. 

```{r}
plt_df %>%
  mutate(facet_title=glue("Coronal slice at:\nworld coord: {coronal_slice_world}\nvoxel coord: {coronal_slice_voxel}")) %>%
  mutate(facet_title=fct_reorder(facet_title, rev(coronal_slice_voxel))) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=intensity), interpolate = T) +
  facet_wrap(~ facet_title, nrow=2) + 
  coord_fixed() + 
  scale_fill_gradient(low="black", high = "white", limits=c(min(arr), 0.75*max_intensity), oob=squish) +
  xlab("Right-left (world) coordinate") +
  ylab("Superior-inferior (world) coordinate") +
  labs(title="Coronal brain slices", 
       subtitle=glue("{length(slices)} evenly spaced slices")) +
  theme_bw() +
  theme(panel.grid = element_blank())
```

## Step 7: Add a slice locator

We'll take contours from a sagittal slice and mark the positions of the coronal slice. 

```{r}
# Get the background file as an array
# Previously, we modified values in this array (setting non-brain values to white)
# Use the unmodified array for contours (so load this in again)
arr_original <- mincArray(mincGetVolume(background_file))

# Pick a sagittal slice for the contour
sagittal_slice <- arr_original[110, , ]

# Pick the intensities at which to determine the contours
contour_levs <- c(20, 50, 80, 100, 120)

# Use grDevices::contourLines() to get the contours
# Returns a list of contour paths
sagittal_contours <- contourLines(x = seq(from=starts[2], to=ends[2], length.out = dim(arr_original)[2]), 
                                  y = seq(from=starts[3], to=ends[3], length.out = dim(arr_original)[3]), 
                                  z = sagittal_slice, 
                                  levels = contour_levs)

# Contour data comes as a list. Put it into a data frame and index each contour path with a new column (obj)
sagittal_contours_df <- 1:length(sagittal_contours) %>%
  map(function(i) {as_tibble(sagittal_contours[[i]]) %>%
      mutate(obj=i)}) %>%
  bind_rows() 

# Many contour paths exist. Remove the small ones by filtering out paths with few points
sagittal_contours_df <- sagittal_contours_df %>%
  group_by(obj) %>%
  filter(length(level) >= 200) %>%
  ungroup()

# Plot contours
sagittal_contours_df %>%
  ggplot(aes(x=x, y=y, group=obj)) +
  geom_path(aes(color=level)) +
  geom_vline(aes(xintercept=world), data=slices_voxel, color='red', lty=1, size=1.5) +
  coord_fixed() +
  scale_color_gradient(name="Intensity level", low = 'grey80', high='grey20') +
  scale_x_continuous(breaks = seq(from=-8, to=5, by=1)) +
  scale_y_continuous(breaks = seq(from=-2, to=5, by=1)) +
  xlab("Anterior-posterior (world) coordinate") +
  ylab("Superior-inferior (world) coordinate") +
  theme_bw()

```


## Step 8:  Combine slice locator with plot and prettify

Finally, put everything together (stylized slice locator and slice plot) using functions from the `patchwork` library. 

```{r}
# Plot the locator (slightly modified from before)
plt_locator <- sagittal_contours_df %>%
  ggplot(aes(x=x, y=y)) +
  geom_path(aes(color=level, group=obj)) +
  geom_segment(aes(x=world, xend=world, y=-2, yend=5), data=slices_voxel, color='red', lty=1, size=0.8) +
  geom_label(aes(x=world, y=-2.5, label=index), 
            hjust=0.5,
            size=3,
            color='red',
            data=slices_voxel) +
  coord_fixed() +
  scale_color_gradient(name="Intensity level", low = 'grey80', high='grey20', guide=F) +
  scale_x_reverse(breaks = c()) +
  scale_y_continuous(breaks = c(), limits=c(-4, 5)) +
  theme_bw() + 
  theme(panel.border = element_blank(),
        panel.grid = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank())

# Plot the slices (slightly modified from before)
plt_slices <- plt_df %>%
  mutate(facet_annotation=glue("AP-coordinate:\n{coronal_slice_world} mm")) %>%
  mutate(facet_annotation=fct_reorder(facet_annotation, rev(coronal_slice_voxel))) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=intensity), interpolate = T) +
  geom_text(aes(x=5, y=6, label=facet_annotation), 
            hjust=1,
            size=3,
            data=. %>% 
              group_by(coronal_slice_voxel) %>%
              slice(1) %>%
              ungroup()) +
  geom_label(aes(x=-4, y=6, label=coronal_slice_index), 
            hjust=0.5,
            size=5,
            color='red',
            data=. %>% 
              group_by(coronal_slice_index) %>%
              slice(1) %>%
              ungroup()) +
  facet_wrap(~ facet_annotation, nrow=2) + 
  coord_fixed() + 
  scale_x_continuous(breaks = seq(from=-5, to=5, by=1)) +
  scale_y_continuous(breaks = seq(from=-2, to=6, by=1), limits = c(-3, 7)) +
  scale_fill_gradient(low="black", high = "white", limits=c(min(arr), 0.75*max_intensity), guide=F, oob=squish) +
  xlab("Right-left (world) coordinate") +
  ylab("Superior-inferior (world) coordinate") +
  labs(title=glue("Coronal brain slices at {length(slices)} evenly spaced coordinates")) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        panel.border = element_blank(),
        strip.text = element_blank())

# Define the combined plot layout for the above two plots
plt_layout <- "
AAAAAAAAAAAAAAAAAABBBBBBBB
AAAAAAAAAAAAAAAAAABBBBBBBB
AAAAAAAAAAAAAAAAAABBBBBBBB
AAAAAAAAAAAAAAAAAA########
AAAAAAAAAAAAAAAAAA########
"

# Combine plots
plt_output <- plt_slices + plt_locator + plot_layout(design = plt_layout)

# Print the final plot
print(plt_output)


```

\newpage
# Comparison of volume differences and brain connectivity

## Goal

Plot volume differences in a Serotonin-transporter knockout mouse model, along with viral tracing data for projections that emanate from the dorsal raphe nuclus (the serotonin center of the brain). 

<div class="alert alert-info">
  <strong>Note:</strong> The previous example showed how to work with world coordinates. The advantage of working with world coordinates is that the plotting does not depend on image resolution, i.e. you can plot images of different resolutions within the same coordinate system and the data should be aligned (if the images themselves are aligned). 
  This time, we will work with voxel coordinates for simplicity. The code presented can be easily modified to work with world coordinates, as done before.
</div>

## Ingredients

Libraries:

```{r}
library(tidyverse)
library(glue)
library(RMINC)
library(ggnewscale)
library(scales)
library(patchwork)
```

Input arguments:

* `background_file`, a path to a .mnc file, corresponding to the background anatomy
* `mask_file`, a path to a .mnc file, corresponding to a brain mask
* `stats_file`, a path to a .mnc file, corresponding to the statistics (the output of `mincLm()`/`mincLmer()`/`mincAnova()` etc..., written out by `mincWriteVolume()`)
* `qvalue_file`, a path to a .mnc file, corresponding to the pvalues adjusted for multiple comparisons ("q-values", the output of `mincFDR()` on the model object)
* `tracer_file`, a path to a .mnc file, corresponding to the viral tracing data of interest
* `slices`, a vector of indices, corresponding to slices in voxel coordinates

```{r}
background_file <- "resources/SERT/SERT_ALL_study_template.mnc"
mask_file <- "resources/SERT/SERT_ALL_study_mask.mnc"
stats_file <- "resources/SERT/SERT_ALL_anova_Genotype_Fstat.mnc"
qvalue_file <- "resources/SERT/SERT_ALL_anova_Genotype_qvalue.mnc"
tracer_file <- "resources/SERT/ABI_projection_density_maxmerged_DR_Slc6a4-Cre_on_SERT_ALL.mnc"
slices <- seq(from=230, to=70, by=-20)
```

## Step 1: Load the data

Load the data.

```{r}
# Read in the MINC volumes
bgvol <- mincGetVolume(background_file)
maskvol <- mincGetVolume(mask_file)
statsvol <- mincGetVolume(stats_file)
qvalvol <- mincGetVolume(qvalue_file)
tracervol <- mincGetVolume(tracer_file)

# Convert to arrays
bgarray <- mincArray(bgvol)
maskarray <- mincArray(maskvol)
statsarray <- mincArray(statsvol)
qvalarray <- mincArray(qvalvol)
tracerarray <- mincArray(tracervol)
```

## Step 2: Generate the contours

Plotting the background anatomy, stats, q-values, and tracer on the same plot can get quite busy. So let's instead use contours for the background, and generate them here.

```{r}
# Pick a sagittal slice for the contour
contours_df <- slices %>%
  map_dfr(
    function(slice_coordinate) {
      
      # Get 2D coronal slices
      coronal_slice_bg <- bgarray[, slice_coordinate + 1, ]
      coronal_slice_mask <- maskarray[, slice_coordinate + 1, ]
      
      # Remove background from non-brain regions, so that there are no contours here
      coronal_slice <- coronal_slice_bg
      coronal_slice[coronal_slice_mask < 0.5] <- 0
      
      # Pick the intensities at which to determine the contours
      contour_levs <- c(200, 500, 1000, 1200, 1250, 1300, 1400, 1450, 1500, 1550, 1600, 1700, 1800)
      
      # Use grDevices::contourLines() to get the contours
      # Returns a list of contour paths
      coronal_slice_contours <- contourLines(x = seq(from=1, to=dim(coronal_slice)[1], length.out = dim(coronal_slice)[1]), 
                                             y = seq(from=1, to=dim(coronal_slice)[2], length.out = dim(coronal_slice)[2]), 
                                             z = coronal_slice, 
                                             levels = contour_levs)
      
      # Contour data comes as a list. Put it into a data frame, index each contour path with a new column (obj), and also index the slice
      coronal_contours_df <- 1:length(coronal_slice_contours) %>%
        map(function(i) {as_tibble(coronal_slice_contours[[i]]) %>%
            mutate(obj=i,
                   slice=slice_coordinate)}) %>%
        bind_rows() 
      
      # Return dataframe
      return(coronal_contours_df)
    }
  )

# Many contour paths exist. Remove the small ones by filtering out paths with few points
contours_df <- contours_df %>%
  group_by(obj, slice) %>%
  filter(length(level) >= 150) %>%
  ungroup()

# Examine the contours
contours_df %>%
  mutate(slice=factor(slice, levels = slices)) %>%
  ggplot(aes(x=x, y=y, group=obj)) +
  geom_path(size=0.2) +
  facet_wrap(~ slice, nrow=3) +
  coord_fixed() +
  scale_color_gradient(name="Intensity level", low = 'grey80', high='grey20') +
  xlab("Right-left (voxel) coordinate") +
  ylab("Superior-inferior (voxel) coordinate") +
  theme_bw()
```

## Step 3: Add the tracer layer

Add the tracer overlay. Also, we'll use `theme_void()` here.

```{r}
# Get the tracer slices into a long form data frame
tracer_df <- slices %>%
  map_dfr(
    function(slice_coordinate) {
      
      # Get 2D coronal slice for tracer projection density
      tracer_slice <- tracerarray[, slice_coordinate + 1, ]
      
      # Label row and column names for this array according to the world coordinates
      rownames(tracer_slice) <- as.character(seq(from=1, to=dim(tracer_slice)[1], length.out = dim(tracer_slice)[1]))
      colnames(tracer_slice) <- as.character(seq(from=1, to=dim(tracer_slice)[2], length.out = dim(tracer_slice)[2]))
      
      # Wrangle data to long form
      out <- tracer_slice %>%
        as_tibble %>%
        mutate(x=rownames(tracer_slice)) %>%
        gather(key = "y", value="projection_density", -one_of("x")) %>%
        mutate(x=as.numeric(x), 
               y=as.numeric(y),
               slice=slice_coordinate)
      
      # Return 
      return(out)
    }
  ) %>%
  mutate(slice=factor(slice, levels = slices)) 

# Plot
contours_df %>%
  mutate(slice=factor(slice, levels = slices)) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=projection_density), 
              interpolate = T,
              data=tracer_df) +
  geom_path(aes(group=obj), size=0.2) +
  facet_wrap(~ slice, nrow=3) +
  coord_fixed() +
  scale_fill_gradient(name="Tracer projection\ndensity",
                      low = "darkgreen", high="green", na.value = "transparent",
                      limits=c(0.1, 0.25)) +
  xlab("Right-left (voxel) coordinate") +
  ylab("Superior-inferior (voxel) coordinate") +
  theme_void()
```

## Step 4: Read in the stats data

Read in the statistics (stats map and FDR-adjusted p-value map) and coerce to long format data frame.

```{r}
# Get stats map into data frame
stats_df <- slices %>%
  map_dfr(
    function(slice_coordinate) {
      
      # Get 2D coronal slice for tracer projection density
      stats_slice <- statsarray[, slice_coordinate + 1, ]
      
      # Label row and column names for this array according to the world coordinates
      rownames(stats_slice) <- as.character(seq(from=1, to=dim(stats_slice)[1], length.out = dim(stats_slice)[1]))
      colnames(stats_slice) <- as.character(seq(from=1, to=dim(stats_slice)[2], length.out = dim(stats_slice)[2]))
      
      # Wrangle data to long form
      out <- stats_slice %>%
        as_tibble %>%
        mutate(x=rownames(stats_slice)) %>%
        gather(key = "y", value="statistic", -one_of("x")) %>%
        mutate(x=as.numeric(x), 
               y=as.numeric(y),
               slice=slice_coordinate)
      
      # Return 
      return(out)
    }
  ) %>%
  mutate(slice=factor(slice, levels = slices)) 

# Get (adjusted) p-value map into data frame
qval_df <- slices %>%
  map_dfr(
    function(slice_coordinate) {
      
      # Get 2D coronal slice for tracer projection density
      qval_slice <- qvalarray[, slice_coordinate + 1, ]
      
      # Label row and column names for this array according to the world coordinates
      rownames(qval_slice) <- as.character(seq(from=1, to=dim(qval_slice)[1], length.out = dim(qval_slice)[1]))
      colnames(qval_slice) <- as.character(seq(from=1, to=dim(qval_slice)[2], length.out = dim(qval_slice)[2]))
      
      # Wrangle data to long form
      out <- qval_slice %>%
        as_tibble %>%
        mutate(x=rownames(qval_slice)) %>%
        gather(key = "y", value="qvalue", -one_of("x")) %>%
        mutate(x=as.numeric(x), 
               y=as.numeric(y),
               slice=slice_coordinate)
      
      # Return 
      return(out)
    }
  ) %>%
  mutate(slice=factor(slice, levels = slices)) 

```

## Step 5: Combine stats and p-value data into a single data frame

Given the relationship between the statistics and p-values, let's combine these two data sources into a single layer, plotting the statistics with a transparency related to the adjusted p-values.

```{r}
# Join data frames
stats_qval_df <- stats_df %>%
  left_join(qval_df, by=c("x", "y", "slice"))

# Add a column that maps qvalue to transparency 
stats_qval_df <- stats_qval_df %>%
  mutate(qalpha=case_when(
    qvalue > 0.2 ~ 0,
    TRUE ~ 0.8*(1 - 5*qvalue))
  )
```


## Step 6: Add the stats layer

Next, add the statistics layer. Note that by default, you cannot have multiple different colour scales for the same aesthetic type. We're using `ggnewscale` to allow for a different colour scales between the tracer and stats layers. 

```{r}
contours_df %>%
  mutate(slice=factor(slice, levels = slices)) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=projection_density), 
              interpolate = T,
              data=tracer_df) +
  scale_fill_gradient(name="Tracer projection\ndensity",
                      low = "darkgreen", high="green", na.value = "transparent",
                      limits=c(0.1, 0.25)) +
  new_scale_fill() +
  geom_raster(aes(fill=statistic,
                  alpha=qalpha), 
              interpolate = T,
              data=stats_qval_df) +
  scale_fill_gradient(name="F-statistic",
                      low = "white", high="red", na.value = "transparent") +
  scale_alpha_continuous(guide=F, limits = c(0,1)) +
  geom_path(aes(group=obj), size=0.2) +
  facet_wrap(~ slice, nrow=3) +
  coord_fixed() +
  xlab("Right-left (voxel) coordinate") +
  ylab("Superior-inferior (voxel) coordinate") +
  theme_void()
```

## Step 7: Add a slice locator

We'll take contours from a sagittal slice and mark the positions of the coronal slice. 

```{r}

# Pick a sagittal slice for the contour
sagittal_slice_bg <- bgarray[110, , ]
sagittal_slice_mask <- maskarray[110, , ]

# Remove background from non-brain regions, so that there are no contours here
sagittal_slice <- sagittal_slice_bg
sagittal_slice[sagittal_slice_mask < 0.5] <- 0

# Pick the intensities at which to determine the contours
contour_levs <- c(200, 500, 1000, 1200, 1250, 1300, 1400, 1450, 1500, 1550, 1600, 1700, 1800)

# Use grDevices::contourLines() to get the contours
# Returns a list of contour paths
sagittal_contours <- contourLines(x = seq(from=1, to=dim(sagittal_slice)[1], length.out = dim(sagittal_slice)[1]), 
                                  y = seq(from=1, to=dim(sagittal_slice)[2], length.out = dim(sagittal_slice)[2]), 
                                  z = sagittal_slice, 
                                  levels = contour_levs)

# Contour data comes as a list. Put it into a data frame and index each contour path with a new column (obj)
sagittal_contours_df <- 1:length(sagittal_contours) %>%
  map(function(i) {as_tibble(sagittal_contours[[i]]) %>%
      mutate(obj=i)}) %>%
  bind_rows() 

# Many contour paths exist. Remove the small ones by filtering out paths with few points
sagittal_contours_df <- sagittal_contours_df %>%
  group_by(obj) %>%
  filter(length(level) >= 150) %>%
  ungroup()

# Plot contours
sagittal_contours_df %>%
  ggplot(aes(x=x, y=y, group=obj)) +
  geom_path(aes(color=level)) +
  geom_vline(aes(xintercept=world), data=tibble(world=slices), color='red', lty=1, size=1.5) +
  coord_fixed() +
  scale_color_gradient(name="Intensity level", low = 'grey80', high='grey20', guide=F) +
  theme_void()

```

## Step 8: Combine slice locator with plot

Combine the above plots (locator and slice series), with additional code to add labels. 

```{r}
# Create a dataframe that maps slice index to the label we want
slice_label_df <- tibble(slice=slices, slice_label=LETTERS[1:length(slices)]) %>%
  mutate(slice=factor(slice, levels=slices),
         slice_label=factor(slice_label, levels=LETTERS[1:length(slices)]))

# Plot the locator (slightly modified from before)
plt_locator <- sagittal_contours_df %>%
  ggplot(aes(x=x, y=y)) +
  geom_path(aes(color=level, group=obj)) +
  geom_vline(aes(xintercept=as.integer(as.character(slice))), data=slice_label_df, color='red', lty=1, size=0.5) +
  geom_label(aes(x=as.integer(as.character(slice)), y=20, label=slice_label), 
             hjust=0.5,
             size=2.5,
             color='red',
             data=slice_label_df) +
  coord_fixed() +
  scale_color_gradient(name="Intensity level", low = 'grey80', high='grey20', guide=F) +
  scale_x_reverse() +
  theme_void()

# Plot the slices (slightly modified from before)
plt_slices <- contours_df %>%
  mutate(slice=factor(slice, levels = slices)) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=projection_density), 
              interpolate = T,
              data=tracer_df) +
  scale_fill_gradient(name="Tracer projection\ndensity",
                      low = "darkgreen", high="green", na.value = "transparent",
                      limits=c(0.1, 0.25)) +
  new_scale_fill() +
  geom_raster(aes(fill=statistic,
                  alpha=qalpha), 
              interpolate = T,
              data=stats_qval_df) +
  scale_fill_gradient(name="F-statistic",
                      low = "white", high="red", na.value = "transparent") +
  scale_alpha_continuous(guide=F, limits = c(0,1)) +
  geom_path(aes(group=obj), size=0.2) +
  geom_label(aes(x=20, y=150, label=slice_label), 
             hjust=0.5,
             size=5,
             color='red',
             data=slice_label_df) +
  facet_wrap(~ slice, nrow=3) +
  coord_fixed() +
  theme_void() +
  theme(strip.text = element_blank())

# Define the combined plot layout for the above two plots
plt_layout <- "
AAAAAAAAAAAAAAAAAABBBBBBBB
AAAAAAAAAAAAAAAAAABBBBBBBB
AAAAAAAAAAAAAAAAAABBBBBBBB
AAAAAAAAAAAAAAAAAA########
AAAAAAAAAAAAAAAAAA########
"

# Combine plots
plt_output <- plt_slices + plt_locator + plot_layout(design = plt_layout)

# Print the final plot
print(plt_output)


```

## Alternative plot: Separately plot data layers

The plot above is still a bit busy, and the transparency (which allows the tracer data to be seen) also might make determining the statistic difficult. So let's separate the data layers and plot them on their own template, similar to the "traditional" plots generated by `mincPlotSliceSeries()` / `sliceSeries()`.

```{r}
# Generate data frame for plotting the background
# Map over each coronal slice
bg_df <- seq_along(slices) %>%
  map_dfr(
    function(i) {
      
      # Get slice
      s <- slices[i]
      
      # Get 2D array corresponding to background intensities
      coronal_bg_array <- bgarray[,s + 1,]
      coronal_mask_array <- maskarray[,s + 1,]
      
      # Fix background 
      coronal_slice_array <- coronal_bg_array
      coronal_slice_array[coronal_mask_array < 0.5] <- 1600
      
      # Label row and column names for this array according to the world coordinates
      rownames(coronal_slice_array) <- as.character(seq(from=1, to=dim(coronal_slice_array)[1], length.out = dim(coronal_slice_array)[1]))
      colnames(coronal_slice_array) <- as.character(seq(from=1, to=dim(coronal_slice_array)[2], length.out = dim(coronal_slice_array)[2]))
      
      # Wrangle data to long form
      out <- coronal_slice_array %>%
        as_tibble %>%
        mutate(x=rownames(coronal_slice_array)) %>%
        gather(key = "y", value="intensity", -one_of("x")) %>%
        mutate(x=as.numeric(x), 
               y=as.numeric(y),
               slice=s)
      
      # Return 
      return(out)
    }
  )

# Plot the stats map
plt_stats <- bg_df %>%
  mutate(slice=factor(slice, levels = slices)) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=intensity), interpolate = T) +
  scale_fill_gradient(low="black", high = "white", limits=c(600, 1600) , oob=squish, guide=F) +
  new_scale_fill() +
  geom_raster(aes(fill=statistic, alpha=alpha),
              interpolate = T,
              data=stats_qval_df %>%
                mutate(alpha=case_when(statistic < 10 ~ 0,
                                       TRUE ~ 0.75))) +
  scale_fill_gradient(name="F-statistic",
                      low = "white", high="red", na.value = "transparent",
                      limits=c(10,50),
                      oob=squish) +
  scale_alpha_continuous(range=c(0,1), guide=F) +
  geom_label(aes(x=20, y=150, label=slice_label), 
             hjust=0.5,
             size=3,
             color='red',
             data=slice_label_df) +
  facet_wrap(~ slice, nrow=3) + 
  coord_fixed() +
  labs(title="F-statistic associated with Slc6a4 genotype") +
  theme_void() +
  theme(strip.text = element_blank())

# Plot the tracer map
plt_tracer <- bg_df %>%
  mutate(slice=factor(slice, levels = slices)) %>%
  ggplot(aes(x=x, y=y)) +
  geom_raster(aes(fill=intensity), interpolate = T) +
  scale_fill_gradient(low="black", high = "white", limits=c(600, 1600) , oob=squish, guide=F) +
  new_scale_fill() +
  geom_raster(aes(fill=projection_density, alpha=alpha), 
              interpolate = T,
              data=tracer_df %>%
                mutate(alpha=case_when(projection_density < 0.1 ~ 0,
                                       TRUE ~ 0.75))) +
  scale_fill_gradient(name="Tracer projection\ndensity",
                      low = "darkgreen", high="green", na.value = "transparent",
                      limits=c(0.1, 0.25),
                      oob=squish) +
  scale_alpha_continuous(range=c(0,1), guide=F) +
  geom_label(aes(x=20, y=150, label=slice_label), 
             hjust=0.5,
             size=3,
             color='red',
             data=slice_label_df) +
  facet_wrap(~ slice, nrow=3) +
  coord_fixed() +
  labs(title="Slc6a4 projections arising from the\nDorsal raphe nucleus") +
  theme_void() +
  theme(strip.text = element_blank())

# Define the combined plot layout for the above two plots
plt_layout <- "
CCCAAAA
CCCAAAA
###AAAA
###AAAA
###BBBB
###BBBB
###BBBB
###BBBB
"

# Combine plots
plt_output <- plt_stats + plt_tracer + plt_locator + plot_layout(design = plt_layout)

# Print the final plot
print(plt_output)
```


<!-- \newpage -->
<!-- # Comparisons of brain atlases -->

<!-- \newpage -->
<!-- # Deformations of an image registration -->

<!-- \newpage -->
<!-- # Visualizing image registration quality -->

<!-- \newpage -->
<!-- # A movie through the brain -->

<!-- \newpage -->
<!-- # A movie of brain registration -->