---
# title: "SMART:\nSemi-manual alignment to reference templates"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# SMART: Semi-manual alignment to reference templates
***
Author(s):
Maintainer: Michelle Jin
- michelle.jin@nih.gov
- <span style="color:#337ab7">@</span>[Michelle_Jin1](https://twitter.com/michelle_jin1)
Contributor: Joseph Nguyen
- jdnguyen@email.wm.edu
### Introduction:
This package builds a pipeline interfacing with the
[wholebrain](https://github.com/tractatus/wholebrain) package ^[Daniel Furth -
@wholebrainsuite (2015). wholebrain: Functions to segment and register cells
from microscopic image files with R. R package version 0.1.1.] to process
whole brain imaging datasets. Datasets of a handful of single slice images
across the brain can also be fed into this pipeline. For brevity, I'll refer
to these datasets as partial brains.
There are five main advancements that SMART builds on top of the wholebrain
base package:
1. Guides the user through setting up analysis parameters and feeds the imaging
data through registration, segmentation, and forward warp process. It automates
looping through these analyses with simple function calls so the user does not
have to set up this looping.
2. Provides a user-friendly interface whenever possible. This includes a console
interface during the registration process, that allows for easy adding,
changing, and removing correspondence points or reverting back to a previous
correspondence change.
3. Helps to solve non-linear relationships between atlas AP coordinates to
actual distances in the anterior-posterior axis in real imaging datasets. This
is critical when dealing with tissue clearing methods that cause non-uniform
tissue morphing across the AP-axis.
4. Organizes and stores the output of wholebrain analysis in automatically
generated subdirectories. This standardized organization will be
beneficial when sharing data output from wholebrain.
5. Provides more ways to parse whole brain datasets and visualize data across
ROIs.
*Disclaimer. This package is intended for coronal plane datasets. Preprocessing
must be done to align the dataset coronally before it should be analyzed with
the pipeline.*
**Below is an illustratratration of the pipeline schematic:**
***
###Installation of SMART:
To use SMART, wholebrain must first be installed. The instructions to install
wholebrain can be found [here](http://www.wholebrainsoftware.org/cms/install).
Please note, if you are using a windows machine, manual installation of opencv
is no longer necessary. I highly recommend reading through some of the documentation
for wholebrain before using this package.
SMART is easily installed from github using the devtools package:
```{r, eval = FALSE}
# Load devtools
library(devtools)
# Install SMART
devtools::install_github("mjin1812/SMART")
# Load SMART
library(SMART)
# Pull up package descriptions and list of functions
?SMART
```
***
### Pipeline breakdown:
As illustrated in the schematic and stated in the package descriptions, not all
parts of this pipeline will be used when analyzing a partial brain. Sections/functions exclusively meant to be used with whole brain datasets will be marked with **(W)**. Additionally, some package functions may be useful for analysis outside the pipeline we present here.
For processing a wholebrain dataset, the pipeline is split into 5 sections. These sections are listed below along with the functions that belong to each section:
**Part 1: Setup pipeline**
+ [`setup_pl()`](#id1) User friendly way to setup parameters for whole or partial brain pipeline analysis.
+ [`im_sort()`](#id2) A function to sort image paths for imaging datasets.
+ [`get_savepaths()`](#id3) Generate savepaths and save directories.
**Part 2: Alignment (W) **
+ [`choice()`](#id4) **(W)** User friendly choice game to align internal reference atlas plates.
+ [`brainmorph()`](#id5) **(W)** Generate a brain morph plot.
+ [`interpolate()`](#id6)**(W)** Interpolate all AP and z numbers for atlas plates in a whole brain project
**Part 3: Registration**
+ [`registration2()`](#id7) Console user interface built on top of registration() function from the wholebrain package.
+ [`regi_loop()`](#id8) Automatically loops through the image registration process for the imaging dataset.
**Part 4: Segmentation, duplicate cleanup, & forward warping**
+ [`filter_loop()`](#id9) Loops through reference slices and allows user to check/change filter settings for segmentation or registration.
+ [`seg_loop()`](#id10) Loop through and segments images in the segmentation channel.
+ [`clean_duplicates()`](#id11) **(W)** Duplicate cell count clean up of segmentation output.
+ [`forward_warp()`](#id12) Performs forward warp loop back onto atlas space using segmentation and registration output.
**Part 5. Dataset manipulation and plotting**
+ [`get_rois()`](#id13) Gets a subset of the forward warped dataframe of just regions of interest.
+ [`get_sunburst()`](#id14) Generates a sunburst plot using a forward warped dataset.
+ [`get_tree()`](#id15) Creates a dataframe of hierarchical region cell count data.
+ [`glassbrain2()`](#id16) Generates a 3D plot of cells counts with option of removing glassbrain in the background.
Below is a walkthrough of each of these functions and the pipeline as a whole
using an example whole brain dataset. Feel free to use the function links above
to skip sections!
For convention, if a return value is given by a function, the recommended variable name is indicated in italics in the return section of each function's help page. If this return value is a required argument for another function in the pipeline, the recommended variable name will also be the same name as the argument.
***
### Part 1: Setup pipeline
This sections sets up analysis parameters, sorted image paths, and generates
savepaths and directories for the rest of the pipeline.
<a id="id1"></a>
**Step 1.**
**`setup_pl()()`** This function asks the user for setup information. Based on input from the user, the function returns a list of parameters for either a whole or partial brain analysis.
```{r, eval = FALSE}
# Scroll to the details section of the help page to see the setup parameters
?setup_pl()
# Run setup_pl, enter parameter information to the console,
# and store the output list into a variable named setup
setup <- setup_pl()
# Check or modify your setup parameters
setup <- setup_pl(setup)
# Note: Whenever a different user works on analyzing the same dataset,
# run the above command to change user initials. This will keep track of
# who did what.
```
> Tips: When providing folder paths, do not put quotes around the path to the console input.
For convention, sequential image numbers are z image numbers, even for a partial brain analysis. Z image number should start at 1 and end at the number of images. Image files should start indexing at 1 in the filenames to match image number.
For a whole brain analysis, the first and last atlas plates must be qualitatively aligned with the first and last z image numbers. Note that the `setup$internal_ref_z`, `setup$regi_AP`, `setup$regi_z` parameters are not user modifiable and will be NULL until the`choice()` or `interpolate()` functions are run.
Additionally, there are a few default pipeline parameters for whole brain analysis:
+ Spacing between registrations (mm). DEFAULT: 0.100
+ Segmentation step (integer). DEFAULT: 1.
+ AP coordinates of internal reference planes. DEFAULT: 1.91, 1.10, -.42, -0.93
, -1.94 , -2.95, -3.96.
The above coordinates correspond to PFC, NAc, antHyp, start of HIP, posterior Hyp, VTA,
and the PAG, respectively. These coordinates will be the atlas plates used to
"calibrate" to internal z images. They'll be used to interpolate and match the z
images of remaining atlas plates. They were chosen because:
1) Based on our experience, they contain easy internal region landmarks that can be consistently
identified by different users
2) They are somewhat evenly spaced throughout the brain.
However, these coordinates are user modifiable to account for user preference
and varying AP ranges of imaging datasets. If you are using this pipeline for
the first time, I recommend you take the default values.
The console code below shows the setup list I am using:
```{r, eval = FALSE}
# Check setup list
> setup <- setup_pl(setup)
$`anim_ID`
[1] "1602"
$user_init
[1] "MJ"
$regi_channel
[1] "D:/SMART_pilot_data/R15_1602_SG_Rein_NoTest_coronal_C02"
$seg_channel
[1] "D:/SMART_pilot_data/R15_1602_SG_Rein_NoTest_coronal_C01"
$output
[1] "D:/SMART_pilot_data/output"
$z_space
[1] 0.0025
$regi_step
[1] 0.1
$seg_step
[1] 1
$first_AP
[1] 2.917519
$first_z
[1] 1
$last_AP
[1] -4.971794
$last_z
[1] 2800
$internal_ref_AP
[1] 1.9060687 1.0969084 -0.4202672 -0.9259924 -1.9374427 -2.9488931 -3.9603435
$internal_ref_z
NULL
$regi_AP
NULL
$regi_z
NULL
Please review your setup information above:
Do you want to change any settings: Y/N?
>
```
<a id="id2"></a>
**Step 2.**
**`im_sort()`** User friendly way to sort images from the registration channel. Asks for user input to account for flexible naming conventions for images.
```{r, eval = FALSE}
# Sort images and store them
# There will be user walk through for the 'separator'
# and 'position' information necessary to sort images
image_paths <- im_sort(setup, extension = "tif")
# Note: Setting the position argument will skip the user walkthrough.
# e.g.
image_paths <- im_sort(setup, extension = "tif", position = 9)
# Check image_paths
image_paths
```
<a id="id3"></a>
**Step 3.**
**`get_savepaths()`** Create standardized save directories based on the setup information from `setup_pl()`. Returns a list of savepaths to subfolders in the output folder.
```{r, eval = FALSE}
# Create savepaths and generate data subdirectories automatically
savepaths <- get_savepaths(setup)
# Check the output folder to see the subdirectories created!
# Show savepaths generated
savepaths
# To save the global environment at any time during an analysis session run:
save.image(savepaths$envir_savepath)
# Going to the R_data folder in the output folder and clicking on the
# .RData folder by date will revert the global environment back to a previous session.
# Saving everyday will prevent you from losing more than a day's worth of work.
```
> Tip: rerun the save.image() expression to update the global environment savepath every day.
Below shows my console output for my savepaths variable:
```{r, eval = FALSE}
# Update and show savepaths
> savepaths <- get_savepaths(setup)
> savepaths
$`envir_savepath`
[1] "D:/SMART_pilot_data/output/R_data/Animal_1602_2019-01-07_MJ.RData"
$out_reference_aligned
[1] "D:/SMART_pilot_data/output/Reference_aligned"
$out_RC_brain_morph
[1] "D:/SMART_pilot_data/output/RC_brain_morph"
$out_auto_registration
[1] "D:/SMART_pilot_data/output/Registrations_auto"
$out_registration
[1] "D:/SMART_pilot_data/output/Registrations_manual"
$out_registration_warps
[1] "D:/SMART_pilot_data/output/Registration_warps"
$out_segmentation_warps
[1] "D:/SMART_pilot_data/output/Segmentation_warpimages"
$out_segmentation_schem
[1] "D:/SMART_pilot_data/output/Segmentation_schematics"
>
```
***
### Part 2: Alignment **(W)**
This section is necessary for whole brain datasets only. Users first "calibrate"
the reference atlas plates with their appropriate z image number by running
`choice()` and playing the choice game. Based on the aligned reference plates,
the user can match all remaining atlas plates with an interpolated z image using
`interpolate().`
<a id="id4"></a>
**Step 4.**
**`choice()` (W) ** User friendly way to play the wholebrain choice game to
align internal reference atlas plates. Automatically saves images of chosen
aligned images in the folder specified by `savepaths$out_reference_aligned`.
Below is my initial console output when I run`choice()` with my `setup`, `savepaths`, and
`image_paths` variables.
```{r, eval = FALSE}
# Run the choice game and save the aligned
# z reference number back into the setup list:
> setup <- choice(setup, savepaths, image_paths, filetype = "tif") # save chosen images as tifs
Play the choice game!
Your reference AP is 1.91
Your current image sampling choice_step is 200
Which image matches best?
Enter 1, 2 or 3:
```
**How the choice game works:**
The choice game will cycle through the internal reference atlas plates
(represented below by dotted vertical lines). Three choice windows will popup
for each reference AP coordinate, corresponding to choices 1,2 or 3,
respectively. The user simply compares the corresponding atlas plate with the 3
choice windows and enters the best qualitative match.
<center>
{width=30%}
</center>
The middlemost window (choice 2) is SMART's best guess at the z image best
aligned with the atlas coordinate based on the first and last aligned atlas
plates (represented above by solid vertical lines).
<center>
{width=60%}
</center>
> Note: The default positions the windows pop up on screen can be set by the
`xpos` argument.
The current `choice_step` indicates how many z images the
anteriormost (choice 1) or posteriormost (choice 3) images are from the middlemost image
(choice 2). Entering 1, 2 or 3 into the console has three outcomes:
1. **Choice 1** The anteriormost image becomes the new choice 2 on the next
choice option.
2. **Choice 2** On the next choice option, there is a smaller `choice_step` on
the left and right images. Choosing this option will progressively "zoom in" on
the `choice_step` options until the steps can't get smaller. The choice game
will then move on to the next atlas coordinate. By default, the smallest step is 10.
3. **Choice 3** The posteriormost image becomes the new choice 2 on the next
choice option.
> Note: the `choice_step` progression is a user modifiable argument.
<center>
{width=30%}
</center>
Once all the internal references atlas coordinates have been aligned, the
`choice()` function will cycle through the midpoint AP coordinates of the
aligned references as a quality check (represented above with black arrows).
After comparing midpoints, any unsatisfactory midpoints become another internal
reference point. The choice game is replayed for those midpoints.
<center>
{width=35%}
</center>
Below are my internal reference AP coordinates and z numbers after
running`choice()` with my setup:
```{r, eval = FALSE}
# Check internal reference AP numbers
> setup$internal_ref_AP
[1] 1.9060687 1.0969084 -0.4202672 -0.9259924 -1.9374427 -2.4431679 -2.9488931 -3.4546183 -3.9603435
[10] -4.4660687
# Check matched image z numbers
> setup$internal_ref_z
[1] 350 657 1305 1425 1813 1978 2202 2422 2541 2681
```
<a id="id5"></a>
**`brainmorph()` (W) ** Generate a brain morph plot showing areas of relative
expansion or compression along the AP axis compared to the Allen Mouse Brain
Atlas AP coodinates (normalized to 1). The reference atlas points used to
generate the morph plot are plotted in red. Setting `saveplot = TRUE` will
save the brain morph into the data folder designated by
`savepaths$out_RC_brain_morph`
```{r, eval = FALSE}
# Generate brainmorph.
brainmorph(setup, savepaths, saveplot = FALSE)
```
<center>
{width=40%}
</center>
<a id="id6"></a>
**Step 5.**
**`interpolate()` (W) ** This function interpolates all corresponding z and AP values for atlas plates that are not reference plates.
```{r, eval = FALSE}
# Interpolate all remaining atlas plates and their z number
setup <- interpolate(setup)
```
Now checking the setup list will show all matched internal AP plates and z numbers:
```{r, eval = FALSE}
> setup$regi_AP
[1] 2.91751908 2.81637405 2.71522901 2.61408397 2.51293893 2.41179389 2.31064885 2.20950382 2.10835878
[10] 2.00721374 1.90606870 1.80492366 1.70377863 1.60263359 1.50148855 1.40034351 1.29919847 1.19805344
[19] 1.09690840 0.99576336 0.89461832 0.79347328 0.69232824 0.59118321 0.49003817 0.38889313 0.28774809
[28] 0.18660305 0.08545802 -0.01568702 -0.11683206 -0.21797710 -0.31912214 -0.42026718 -0.52141221 -0.62255725
[37] -0.72370229 -0.82484733 -0.92599237 -1.02713740 -1.12828244 -1.22942748 -1.33057252 -1.43171756 -1.53286260
[46] -1.53286260 -1.63400763 -1.73515267 -1.83629771 -1.93744275 -2.03858779 -2.13973282 -2.24087786 -2.34202290
[55] -2.44316794 -2.54431298 -2.64545802 -2.74660305 -2.84774809 -2.94889313 -3.05003817 -3.15118321 -3.25232824
[64] -3.35347328 -3.45461832 -3.55576336 -3.65690840 -3.75805344 -3.85919847 -3.96034351 -4.06148855 -4.16263359
[73] -4.26377863 -4.36492366 -4.46606870 -4.56721374 -4.66835878 -4.76950382 -4.87064885 -4.97179389
> setup$regi_z
[1] 1 36 71 106 141 176 210 245 280 315 350 388 427 465 504 542 580 619 657 700 743
[22] 787 830 873 916 959 1003 1046 1089 1132 1175 1219 1262 1305 1329 1353 1377 1401 1425 1464 1503 1541
[43] 1580 1619 1658 1658 1697 1735 1774 1813 1846 1879 1912 1945 1978 2023 2068 2112 2157 2202 2246 2290 2334
[64] 2378 2422 2446 2470 2493 2517 2541 2569 2597 2625 2653 2681 2705 2729 2752 2776 2800
```
***
### Part 3: Registration
This section automates looping through through and registering all atlas plates
to images from the registration channel. It also allows for easy modification of
registrations.
<a id="id7"></a>
**`registration2()`** Provides a user friendly interface to alter registrations. It
retains the same arguments from the `registration()` function in wholebrain.
This function may be useful for those who already have their own established
analysis pipelines. Check out the function help page for more details.
<center>
{width=70%}
</center>
The schematic above illustrates the user-friendly console interface that allows
for modification of the correspondence plates on a given atlas plate. The
`regi_loop()` function integrates this user-friendly registration interface with
the pipeline.
<a id="id8"></a>
**Step 6.**
**`regi_loop()`** User friendly function to alter atlas to image registration.
When run with defaults, this function will automatically create a list in the
global environment called `regis`. It will loop through the atlas plates
according to the user preferences. The user also has the option to save the
global environment after every manual registration.
```{r, eval = FALSE}
# Example run of automated registration loop
# Create filter for registration
filter<-structure(list(alim = c(50, 50),
threshold.range = c(50000, 60000L),
eccentricity = 999L,
Max = 1000,
Min = 150,
brain.threshold = 270L,
resize = 0.25,
blur = 7L,
downsample = 2))
# Run the autoloop
regi_loop(setup, savepaths, image_paths, autoloop = TRUE, filter = filter)
```
> Note: The function [`filter_loop()`](#id17) allows the user to modify the
filter for registration or segmentation and loop through the references plates
to check filter adjustments in images across the brain.
There are four different loop options:
1. `autoloop = TRUE` runs wholebrain's first pass at registering all atlas plates.
Registrations are stored in the *Registrations_auto* folder. This is useful to
check the quality of registrations across all atlas plates using current filter
settings.
2. `reference = TRUE` loops through reference atlas plates only and displays the
console interface.
3. The `touchup` argument is designed to edit a subset of atlas plates that have
been previously registered
+ `touchup = TRUE` runs a user friendly interface allowing user to enter which
plates they want to fix.
+ `touchup =` Numeric vector of plate numbers (integers) the user wants to
fix.
4. Setting all the arguments above to `FALSE` will loop through all atlas plates
and display the console interface.
>Tip: Save the global environment after every 2-3 registrations. The save time
lengthens significantly as more registrations are stored in the `regis` list.
Once the `regi_loop()` function has been run for the first time, the `regis`
argument must be set to the `regis` list automatically generated in the global
environment.
For example:
```{r, eval = FALSE}
# Loop through and modify only reference atlas plates
regi_loop(setup, savepaths, image_paths, regis = regis, filter = filter, reference = TRUE)
```
Below is an example of a registration before and after manual improvement:
<center>
{width=50%}
{width=50%}
</center>
***
### Part 4: Segmentation, duplicate cleanup, & forward warping
This section loops through all the images of the segmentation channel and
automates the forward looping process. At the end of this section, you will have
a mapped dataset!
<a id="id9"></a>
<a id="id17"></a>
**Step 7.**
**`filter_loop()`** This function takes a registration or segmentation filter
and loops through all images in a partial brain or just the reference images in
either the registration or segmentation channel. The `wholebrain::segment()`
function is used to display the segmentation output. There is a console menu
option for the user to adjust filter parameters if they need.
```{r, eval = FALSE}
# Check and modify filter settings in the segmentaion channel
filter <- filter_loop(setup, image_paths, channel = "seg")
```
<a id="id10"></a>
**Step 8.**
**`seg_loop()`** This function loops through the z images in the segmentation folder path. For a whole brain, every N images is segmented, where N is equal to `setup$seg_step`. For a partial brain, all the images in `setup$regi_z` are segmented.
```{r, eval = FALSE}
# Run segmentation loop and store the output
segs <- seg_loop(setup, image_paths, filter)
```
> Note: If the imaging dataset is large, steps 8, 9, & 10 will be time intensive processes. Processing time will be printed after each of these functions are finished.
<a id="id11"></a>
**Step 9.**
**`clean_duplicates()` (W) ** Duplicate cell count clean up of segmentation output from the `seg_loop()` function. This step is only necessary for large imaging datasets where the z step size is smaller than the width of a cell body.
```{r, eval = FALSE}
# Clean duplicate cell counts
segs <- clean_duplicates(setup, segs, z_thresh = 10, compare_depth = 200)
# z_thresh - Threshold (in um) in the z-plane for the maximum z distance
# before a cell is counted as another cell.
# compare_depth - Comparision depth in (um). Compare a single image with adjacent
# images up to 200 um posterior.
```
<a id="id12"></a>
**Step 10.**
**`forward_warp()` (W) ** Loop through all segmented cells and perform forward
warp loop back onto atlas space using registration and segmentation data. There
are user options for plotting a schematic plot and saving forward warp images.
The function will automatically interpolate AP coordinate values of z images in
between atlas plates and assign them to the dataset.
```{r, eval = FALSE}
# Perform forward warp
dataset <- forward_warp(setup, savepaths, segs, regis)
```
Below are representative forward warp and schematic images:
<center>
**Forward warp image**
{width=75%}
**Schematic plot**
{width=40%}
</center>
> Note: `forward_warp()` will automatically clean up any cell counts that aren't assigned a region ID.
***
### Part 5: Dataset manipulation and plotting:
Following Step 10, there are no more steps in the pipeline. The remaining
functions are designed for easy manipulation and plotting of the dataset.
<a id="id13"></a>
**`get_rois()`** Allows the user to enter a character vector of Allen Mouse
Brain Atlas abbreviations of regions of interest (ROIs). A subset of the
dataframe from the wholebrain dataset of just the ROIs (and their subregions)
are returned.
`get_rois()` is especially useful in plotting cell count tables of regions of interest:
```{r, eval = FALSE}
# Get dataset of just rois
rois <- get_rois(dataset, rois = c("ACB", "CEA", "BLA", "HIP"))
# Plot region cell counts in a table
quartz()
wholebrain::dot.plot(rois)
```
<center>
{width=70%}
</center>
<a id="id14"></a>
**`get_sunburst()`** Generate a sunburst plot using a forward warped dataset.
```{r, eval = FALSE}
# Plot sunburst plot using the get_sunburst() function
SMART::get_sunburst(dataset)
```
```{r, echo = FALSE, cache = FALSE}
load("C:/Users/mjin1/Documents/Test/packagetest_W/R_data/dataset.RData")
SMART::get_sunburst(dataset)
```
> Note: By setting the `rois` argument, the sunburst will display only ROI data.
If the `parent` argument is set to `FALSE`, the base layer in sunburst will be
the first ROI.
<a id="id15"></a>
**`get_tree()`** Create a dataframe of hierarchical region cell count data.
The code below generates a dataframe of the hierarchy tree for the hypothalamus:
```{r, eval = FALSE}
# Hierarchical dataframe of just the hypothalamus
tree <- get_tree(dataset, rois = "HY")
View(tree)
```
> Note: this dataframe may be useful for generating other heirarchical plots such as [treemaps](https://github.com/cran/treemap).
\newpage
<a id="id16"></a>
**`glassbrain2()` ** A modified version of wholebrain::glassbrain(). New options include:
1. Gives the user an option to turn off the "jitter" in cell placement in original glassbrain function (the jitter gives a more space filled look).
2. Does not show the glass brain display when `high.res` is set to "OFF" (otherwise set to TRUE or FALSE). Setting high.res to TRUE or FALSE will render a high or low resolution atlas boundary
```{r, eval = FALSE}
# Generate rgl object and plot glassbrain without 3D atlas boundaries
glassbrain <- glassbrain2(dataset, high.res = "OFF", jitter = FALSE)
# visualize glass brain
glassbrain
```
```{r, echo = FALSE, cache = FALSE}
htmltools::includeHTML(file.path(getwd(), "/schematics/glassbrain.html"))
```
> Tip: Combine the `get_rois()` function output with `glassbrain2()` to get a glass brain of just ROIs.
matiasandina/SMART
Semi-manual alignment to reference templates (SMART): An open source pipeline in R for whole brain mapping
R