paper_code.R

In adigherman/Rxnat: Queries and Extracts Images from Extensible Neuroimaging Archive Toolkit Public/Private Datasets

# Install Rxnat dependencies for the examplein the paper
# The paper example requires cmake, ITK and FSL, as a system dependencies
# cmake installation guide can be found here: https://cmake.org/install/
# To install FSL please follow this guide: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FslInstallation
# The ITK installation guide can be found here: https://itk.org/download/
#
# The R dependencies can be installed easily using 
# `install.packages()` function 
# install.packages(c('dplyr','scales'))
#
# For the `fslr`, `extrantsr`, `malf.templates` and `WhiteStripe`
# we recomment using the Neuroconductor platform to easily install them
# source("https://neuroconductor.org/neurocLite.R")
# neuro_install('fslr')
# neuro_install('extrantsr')
# neuro_install('malf.templates')
# neuro_install('WhiteStripe')
#
# Instead of installing all the system dependencies and package dependencies 
# from above, one can easily download the neurodocker image from this
# page: https://neuroconductor.org/docker-release
#
# Authentication against the XNAT server can be done in two
# different ways by setting up environment variables in R, or 
# by providing the user credentials in the `xnat_connect()`
# function. The recommended and more secure method is by
# setting up the environment variables in the .Renvion file.
# The .Renviron file is located in the user folder (eg. ~/.Renviron)
# and the variable names for the XNAT connection need to be in the 
# form NITRC_RXNAT_USER and NITRC_RXNAT_PASS (for NITRC connection).
# As an example, we can set NITRC_RXNAT_USER="testuser" and 
# NITRC_RXNAT_PASS="testpassword" in the .Renviron file, start 
# (re-start) R and initiate the NITRC connection with 
# nitrc <- xnat_connect("https://nitrc.org/ir", xnat_name="NITRC")
# For any other XNAT server, please replace the first part of the
# environment variable with the XNAT server name (for example
# HCP will have HCP_RXNAT_USER and HCP_RXNAT_PASS). Please make
# sure that the `xnat_name` parameter in the `xnat_connect` call
# is using the same name as the first part of the environment variables
# (eg. NITRC, HCP).

library(dplyr)
library(fslr)
library(extrantsr)
library(malf.templates)
library(scales)
library(WhiteStripe)

##install.packages('Rxnat')

# Establish XNAT connections
library(Rxnat)
nitrc <- xnat_connect("https://nitrc.org/ir", xnat_name="NITRC")
hcp <- xnat_connect("https://db.humanconnectome.org", xnat_name = "hcp")

# Get list subjects from XNAT
nitrc_subjects <- nitrc$subjects()
hcp_subjects <- hcp$subjects()
all_subjects <- bind_rows(nitrc_subjects, hcp_subjects)
colnames(all_subjects)

# Extract experiment data
nitrc_T1_scan_resources <- nitrc$scans(type="T1")
nitrc_T1_scan_resources <- nitrc_T1_scan_resources %>%
  mutate(source = "nitrc", scantype = "T1")
hcp_T1_scan_resources <- hcp$scans(type="T1w")
hcp_T1_scan_resources <- hcp_T1_scan_resources %>%
  mutate(source = "hcp", scantype = "T1")

# Filter resources and select subjects between 26 and 40 years
T1_resources <- bind_rows(nitrc_T1_scan_resources, hcp_T1_scan_resources)
T1_resources <- left_join(
  T1_resources,
  all_subjects,
  by = c("subject_ID" = "ID")
)
age_26_to_40_group <- T1_resources %>%
  filter(Age>26) %>%
  filter(Age<40)


# Download the first subject T1 weighted image
file_path <- nitrc$download_dir(
  experiment_ID = age_26_to_40_group$experiment_ID[8],
  scan_type = "T1",
  extract = TRUE
)


# Read T1 image
t1 <- readrpi(file_path[1])
ortho2(t1, add.orient = TRUE)

# Remove neck and drop empty dimensions
noneck = remove_neck( file_path,
                      template.file = fslr::mni_fname(brain = TRUE, mm = 1),
                      template.mask = fslr::mni_fname(mm = 1, brain = TRUE, mask = TRUE)
)
red = dropEmptyImageDimensions(noneck)
red <- readrpi(red)


# Inhomegeneity correction
t1_n4 = bias_correct(red,
                     correction = "N4",
                     outfile = tempfile(fileext = ".nii.gz"), retimg = FALSE
)
t1_n4 <- readrpi(t1_n4)

# Perform registration against malf
timgs = mass_images(n_templates = 35)
ss = malf(infile = t1_n4,
          template.images = timgs$images,
          template.structs = timgs$masks,
          keep_images = FALSE)

# Perform skull stripping
proc_outfile <- paste0("T1_Processed.nii.gz")
proc_outfile <- file.path(tempdir(),proc_outfile)
skull_ss <- preprocess_mri_within(
  files = t1_n4,
  outfiles = proc_outfile,
  correction = "N4",
  maskfile = ss,
  correct_after_mask = FALSE)
t1_ss <- readrpi(proc_outfile)

# Perform WhiteStripe intensity normalization
ind = whitestripe(img = t1_ss, type = "T1", stripped = TRUE)$whitestripe.ind
ws_t1 = whitestripe_norm(t1_ss, ind)

# Perform segmentation
ss_tcs = fslr::fast_nobias(ws_t1,
                           verbose = TRUE)

# Generate Panel A
ortho2(red,
       add.orient = TRUE,
       addlegend = TRUE,
       leg.cex = 1.5,
       leg.title = "Panel A",
       leg.col = "red",
       legend = "T1 Image Neck Removed",
       leg.x = 0
)

# Generate Panel B
ortho2(red,
       t1_ss,
       col.y=alpha("red", 0.3),
       add.orient = TRUE,
       addlegend = TRUE,
       leg.cex = 1.5,
       leg.title = "Panel B",
       leg.col = "red",
       legend = "T1 Image Brain Mask",
       leg.x = 0
)

# Generate Panel C
double_ortho(ws_t1,
             ss_tcs,
             add.orient=TRUE,
             addlegend =TRUE,
             leg.cex=1.5,
             leg.title="Panel C",
             leg.col="red",
             legend="Segmentation results",
             leg.x=0
)