Data preparation

Libraries

library(RMINC)
## Loading required package: BatchJobs
## Loading required package: BBmisc
## Sourcing configuration file: '/axiom2/projects/software/arch/linux-xenial-xerus/RMINC/1.5.0.0/BatchJobs/etc/BatchJobs_global_config.R'
## BatchJobs configuration:
##   cluster functions: Interactive
##   mail.from: 
##   mail.to: 
##   mail.start: none
##   mail.done: none
##   mail.error: none
##   default.resources: 
##   debug: FALSE
##   raise.warnings: FALSE
##   staged.queries: TRUE
##   max.concurrent.jobs: Inf
##   fs.timeout: NA
## Sourcing configuration file: '/micehome/chammill/R/x86_64-pc-linux-gnu-library/3.3/RMINC/parallel/sge_BatchJobs.R'
library(ggplot2)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:stats':
## 
##     sigma

Read in bits of data

## Set a data path
data_path <- "/hpf/largeprojects/MICe/chammill/presentations/summer_school2017/longitudinal_analysis/comparing_reg_chain_and_tamarack/"

gf <- read.csv(file.path(data_path, "../tamarack_det.csv"))
load(file.path(data_path, "../anatAbsVols.RData"))

Set anatomy definitions

defsFile <- "/hpf/largeprojects/MICe/tools/atlases/Dorr_2008_Steadman_2013_Ullmann_2013/mappings/Dorr_2008_Steadman_2013_Ullmann_2013_mapping_of_labels.csv"

And make a hierarchical tree

AllenDefs <- "/hpf/largeprojects/MICe/tools/atlases/Allen_Brain/Allen_hierarchy_definitions.json"
hdefs <- makeMICeDefsHierachical(defsFile, AllenDefs)

Add volumes

hanat <- addVolumesToHierarchy(hdefs, anatAbsVols)

Slope type

library(broom)
library(splines)
library(data.tree)
slopeComp <- Clone(hanat)
slopeComp$Do(function(brainROI) {
  gf$volume <- brainROI$volumes
  l1 <- lmer(volume ~ ns(timepoint, 1) + (1|mouse), data=gf, REML=FALSE)
  l2 <- lmer(volume ~ ns(timepoint, 2) + (1|mouse), data=gf, REML=FALSE)
  l3 <- lmer(volume ~ ns(timepoint, 3) + (1|mouse), data=gf, REML=FALSE)
  brainROI$l1v2 <- tidy(anova(l1, l2))$statistic[2]
  brainROI$l2v3 <- tidy(anova(l2, l3))$statistic[2]
})
hanatView(slopeComp, "l2v3", low=10, high=1000)
gf$volume <- FindNode(hanat, "Posterior parietal association areas")$volumes
qplot(timepoint, volume, data=gf) + geom_smooth(method="lm", formula=y~ns(x,2)) +
  geom_smooth(method="lm", formula=y~ns(x,3), col=I("red"))

anatVolume <- mincArray(mincGetVolume("/hpf/largeprojects/MICe/chammill/presentations/summer_school2017/longitudinal_analysis/DSU_atlas/Dorr_2008_Steadman_2013_Ullmann_2013_on_MEMRI_C57BL6_P43_average.mnc"))
labelVolume <- mincArray(mincGetVolume("/hpf/largeprojects/MICe/chammill/presentations/summer_school2017/longitudinal_analysis/DSU_atlas/Dorr_2008_Steadman_2013_Ullmann_2013_on_MEMRI_C57BL6_P43_labels.mnc"))

vol <- hanatToVolume(slopeComp, labelVolume, "l2v3")
mincPlotSliceSeries(anatVolume, vol, low=10, high=100, anatLow=700, anatHigh=1400, begin=50, end=-50)

chivals <- slopeComp$Get("l2v3")
qvals <- p.adjust(pchisq(chivals, 2, lower.tail=F), "BY")
slopeComp$Set(l2v3q = qvals)
slopeComp$Set(l2v3ql10 = -log10(qvals))
hanatView(slopeComp, "l2v3ql10", low=1.3, high=4)

Individual growth patterns

indComp <- Clone(hanat)
indComp$Do(function(x) {
  gf$volume <- x$volumes
  gfwt <- subset(gf, genotype=="WT" & timepoint <= 12)
  l1 <- lmer(volume ~ ns(timepoint, 2) + (1|mouse), data=gfwt, REML=F)
  l2 <- lmer(volume ~ ns(timepoint, 2) + (timepoint|mouse), data=gfwt, REML=F)
  x$indvgroup <- tidy(anova(l1, l2))$statistic[2]
})
hanatView(indComp, "indvgroup", low=3, high=10)
vol <- hanatToVolume(indComp, labelVolume, "indvgroup")
mincPlotSliceSeries(anatVolume, vol, low=2, high=10, anatLow=700, anatHigh=1400, begin=50, end=-50)

indComp2 <- Clone(indComp)
Prune(indComp2, function(x) !startsWith(x$name, "left") & !startsWith(x$name, "right"))
## [1] 262
vol <- hanatToVolume(indComp2, labelVolume, "indvgroup")
mincPlotSliceSeries(anatVolume, vol, low=2, high=10, anatLow=700, anatHigh=1400, begin=50, end=-50)