MBM_test
Nick Wang 2018-08-20
This RMarkdown notebook serves as a test bed for pydpiper's canonical MBM.py run. To run this Rmd file as a script on the command line: Rscript -e "rmarkdown::render('MBM_test.Rmd')"
Download the required data from the web.
wget repo.mouseimaging.ca/repo/Pydpiper_test_files/test-data_20180925.tar.gz
tar -xjf test-data_20180925.tar.gz test-data_20180925Running MBM.py on the command line - this no longer works!!!
Run tail -F MBM_out.txt to track the standard output of the following MBM.py call. This MBM.py call runs a 6 parameter alignment towards an initial model assuming the input files have a random orientation scanned in different coils/spaces (lsq6-large-rotations), then runs a 12 parameter alignment towards a linear consensus average, then a non-linear alignment towards a non-linear consensus average; MAGeT segmentation is done on the lsq6 files to segment the brains.
#Run `tail -F MBM_out.txt` to track the standard output
MBM.py --pipeline-name=MBM_test \
--subject-matter mousebrain \
--num-executors 1000 --time 48:00:00 \
--csv-file test-data_20180925/input.csv \
--lsq6-large-rotations-tmp-dir=/tmp \
--init-model test-data_20180925/basket_mouse_brain_40micron.mnc \
\
--run-maget \
--maget-registration-method minctracc \
--maget-atlas-library test-data_20180925/ex-vivo/ \
--maget-nlin-protocol test-data_20180925/default_nlin_MAGeT_minctracc_prot.csv \
--maget-masking-nlin-protocol test-data_20180925/default_nlin_MAGeT_minctracc_prot.csv \
\
--lsq12-protocol test-data_20180925/Pydpiper_testing_default_lsq12.csv \
>> MBM_out.txtThe canonical MBM.py analysis uses functions from tidyverse, RMINC, and visualization tools in MRIcrotome.
library(tidyverse)## ── Attaching packages ────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.0.0 ✔ purrr 0.2.5
## ✔ tibble 1.4.2 ✔ dplyr 0.7.6
## ✔ tidyr 0.8.1 ✔ stringr 1.3.1
## ✔ readr 1.1.1 ✔ forcats 0.3.0
## ── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(RMINC)
library(grid)
library(MRIcrotome)##
## Attaching package: 'MRIcrotome'
## The following object is masked from 'package:dplyr':
##
## slice
## The following object is masked from 'package:graphics':
##
## legend
#This is a fix for knittr
knitr::opts_chunk$set(cache = TRUE, warning = FALSE, message = FALSE, cache.lazy = FALSE)In total, there are 16 brains imaged twice, 2 weeks apart. The second set of images had volume changes artificially induced in the following regions:
## # A tibble: 5 x 2
## region inducedChange
## <chr> <dbl>
## 1 olfactory bulbs 1.05
## 2 striatum 0.85
## 3 cerebral cortex: occipital lobe 0.93
## 4 dentate gyrus of hippocampus 1.1
## 5 cerebellar cortex 1.07
Load the consensus average and corresponding masks as mincArrays.
#mincArray gives the mincSingleDim appropriate dimensions
consensusVol <- file.path("MBM_test_nlin", "nlin-3.mnc") %>%
mincGetVolume() %>%
mincArray
d <- consensusVol %>% dim()
consensusMaskPath <- file.path('MBM_test_nlin', 'nlin-3_mask.mnc')
consensusMask <- consensusMaskPath %>% mincGetVolumeLoad MBM.py's useful output csv file pointing to the processed images.
gfs <- file.path("analysis.csv") %>%
read_csv() %>%
filter(fwhm == 0.2) %>%
mutate(name = str_c(group, coil, sep="_"),
group = fct_relevel(group, "wt")) %>%
select(-fwhm)Do an anatGetAll call on the MAGeT segmented label files. anatMatrix is a matrix of volumes, one for each region for each brain.
defPath <- "test-data_20180925/ex-vivo/Dorr_2008_mapping_of_labels.csv"
anatMatrix <- anatGetAll(gfs$label_file, method="labels", defs=defPath) %>%
anatCombineStructures(defs=defPath) %>%
unclass() %>%
as_tibble()Find the mean volumes for each region for the two groups.
fractions <- anatMatrix %>%
bind_cols(select(gfs,group)) %>%
gather(key=region, value=size, -group) %>%
group_by(group, region) %>%
summarise(mean_size=mean(size)) %>%
spread(key=group, value=mean_size) %>%
#filter(wt > 1) %>%
mutate(fraction = mut/wt) %>%
arrange(desc(fraction)) %>%
left_join(changes)
fractions## # A tibble: 62 x 5
## region wt mut fraction inducedChange
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 lateral ventricle 4.33 4.80 1.11 NA
## 2 stria terminalis 0.905 0.956 1.06 NA
## 3 cerebellar cortex 44.0 46.4 1.05 1.07
## 4 dentate gyrus of hippocampus 3.76 3.95 1.05 1.1
## 5 stratum granulosum of hippocamp… 0.942 0.984 1.04 NA
## 6 olfactory bulbs 26.0 27.1 1.04 1.05
## 7 arbor vita of cerebellum 10.1 10.5 1.04 NA
## 8 subependymale zone / rhinocele 0.0684 0.0707 1.03 NA
## 9 fimbria 3.29 3.40 1.03 NA
## 10 cerebral aqueduct 0.620 0.635 1.03 NA
## # ... with 52 more rows
Also look at the regions with greatest shrinkage.
top_n(fractions, -10, fraction)## # A tibble: 10 x 5
## region wt mut fraction inducedChange
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 nucleus accumbens 4.15 4.11 0.991 NA
## 2 superior olivary complex 0.764 0.758 0.991 NA
## 3 globus pallidus 3.16 3.13 0.991 NA
## 4 internal capsule 2.88 2.85 0.989 NA
## 5 olfactory tubercle 3.50 3.45 0.985 NA
## 6 corpus callosum 16.0 15.6 0.978 NA
## 7 cerebral cortex: occipital lobe 5.67 5.52 0.973 0.93
## 8 habenular commissure 0.0340 0.0331 0.973 NA
## 9 anterior commissure: pars poste… 0.439 0.426 0.971 NA
## 10 striatum 20.5 18.2 0.891 0.85
Are these changes significant? Call anatLm to fit a linear model on each brain region, and anatFDR to correct for multiple comparisons. Did modified regions yield significant results?
####
anatMatrix <- anatMatrix %>% select(fractions$region)
avs <- anatLm(~group, gfs, anatMatrix)## N: 32 P: 2
## Beginning vertex loop: 62 3
## Done with vertex loop
qavs <- anatFDR(avs)##
## Computing FDR threshold for all columns
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-groupmut
fractions$qvalue_abs <- qavs[,"qvalue-tvalue-groupmut"]
fractions %>%
arrange(qvalue_abs) %>%
filter(!is.na(inducedChange))## # A tibble: 5 x 6
## region wt mut fraction inducedChange qvalue_abs
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 striatum 20.5 18.2 0.891 0.85 0.0000191
## 2 dentate gyrus of hippocam… 3.76 3.95 1.05 1.1 0.00558
## 3 cerebellar cortex 44.0 46.4 1.05 1.07 0.00583
## 4 olfactory bulbs 26.0 27.1 1.04 1.05 0.108
## 5 cerebral cortex: occipita… 5.67 5.52 0.973 0.93 0.748
Repeat the analysis for each region's volume relative to the total brain size by co-varying for brain size in the anatLm call. Do modified regions yield significant results when covarying for total brain size?
gfs$brainVolumes <- anatMatrix %>%
rowSums()
avsrel <- anatLm(~group+brainVolumes, gfs, anatMatrix)## N: 32 P: 3
## Beginning vertex loop: 62 4
## Done with vertex loop
qavsrel <- anatFDR(avsrel)##
## Computing FDR threshold for all columns
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-groupmut
## Computing threshold for tvalue-brainVolumes
fractions$qvalue_rel <- qavsrel[,"qvalue-tvalue-groupmut"]
fractions %>%
arrange(qvalue_rel) %>%
filter(!is.na(inducedChange))## # A tibble: 5 x 7
## region wt mut fraction inducedChange qvalue_abs qvalue_rel
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 striatum 20.5 18.2 0.891 0.85 0.0000191 6.14e-13
## 2 cerebellar cor… 44.0 46.4 1.05 1.07 0.00583 3.33e- 3
## 3 dentate gyrus … 3.76 3.95 1.05 1.1 0.00558 3.33e- 3
## 4 cerebral corte… 5.67 5.52 0.973 0.93 0.748 6.77e- 2
## 5 olfactory bulbs 26.0 27.1 1.04 1.05 0.108 8.18e- 2
For the sake of visualization, let us do voxel-wise analysis.
vs <- mincLm(log_full_det ~ group, gfs, mask=consensusMaskPath)## Method: lm
## Number of volumes: 32
## Volume sizes: 241 478 315
## N: 32 P: 2
## In slice
## 0 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 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## Done
vsFDR <- mincFDR(vs, mask=consensusMaskPath, method="FDR")##
## Computing FDR threshold for all columns
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-groupmut
vsFDR %>% thresholds()## F-statistic tvalue-(Intercept) tvalue-groupmut
## 0.01 25.233605 3.558889 5.023306
## 0.05 15.578418 2.728478 3.946951
## 0.1 12.256784 2.320473 3.500969
## 0.15 10.290461 2.057694 3.207875
## 0.2 8.896389 1.855990 2.982682
vsrel <- mincLm(log_full_det ~ group+brainVolumes, gfs, mask=consensusMaskPath)## Method: lm
## Number of volumes: 32
## Volume sizes: 241 478 315
## N: 32 P: 3
## In slice
## 0 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 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## Done
vsrelFDR <- mincFDR(vsrel, mask=consensusMaskPath, method="FDR")##
## Computing FDR threshold for all columns
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-groupmut
## Computing threshold for tvalue-brainVolumes
vsrelFDR %>% thresholds()## F-statistic tvalue-(Intercept) tvalue-groupmut tvalue-brainVolumes
## 0.01 6.964716 3.150880 4.487269 3.145179
## 0.05 4.250595 2.352565 3.607670 2.349146
## 0.1 3.191952 1.968927 3.207294 1.966168
## 0.15 2.597883 1.725112 2.949397 1.722711
## 0.2 2.187409 1.539876 2.755862 1.537614
sliceSeries(nrow=5, ncol=2, begin=100, end =300) %>%
anatomy(consensusVol, range(consensusVol)[1], range(consensusVol)[2]) %>%
overlay(mincArray(vs, "tvalue-groupmut"), low=vsFDR %>% thresholds() %>% {.["0.05", "tvalue-groupmut"]}, high=10, symmetric = TRUE) %>%
addtitle("Absolute Volume Changes") %>%
contourSliceIndicator(consensusVol, c(700,1400)) %>%
legend("t-statistics") %>%
#######
sliceSeries(nrow=5, ncol=2, begin=100, end = 300) %>%
anatomy() %>% #reuse previous anatomy call's arguments
overlay(mincArray(vsrel, "tvalue-groupmut"), low=vsrelFDR %>% thresholds() %>% {.["0.05", "tvalue-groupmut"]}, high=10, symmetric = TRUE) %>%
addtitle("Relative Volume Changes") %>%
contourSliceIndicator(consensusVol, c(700,1400)) %>%
legend("t-statistics") %>%
draw()