-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathpaper_code.R
170 lines (150 loc) · 5.4 KB
/
paper_code.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
# 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
)