-
Notifications
You must be signed in to change notification settings - Fork 10
Registration chain
Matthijs van Eede edited this page Oct 14, 2016
·
7 revisions
The registration chain is a image registration pipeline perfectly suited to deal with longitudinal data.
TODO:
- add information about what the input is in this example
- explain difference first level and second level (vs. timepoint/subject)
- how do we create a mask
- note about Darren vs RC depending on the input data.
- use Chris's new threshold function to get thresholds instead of [3,3] ...
- add note on how to generate atlases below (and add atlas generation to the registration chain ...)
- make this HOWTO into a Sweave/Rmarkdown document and provide some downloadable analysis files
When your pipeline is finished, a comma separated values (csv) file is generated that can directly be used in R to analyze your data.
...
The first thing you might want to do is analyze the differences at each of the time points in your data (in this example, we'll assume that the first level is within-subject registrations, so the 'timepoint' column of your input CSV really corresponds to the timepoint ...). You can do that as follows:
library(dplyr)
library(RMINC)
# rc is short for registration chain
# load the analysis csv
rc_data_temp <- read.csv("[pipe_name]_analysis_files.csv")
# combine the csv file that was used for the registration chain pipeline and the analysis file that was produced
rc_input_csv <- read.csv("registration_chain_csv_file.csv")
# our data frame does not contain the male/female sex attribute, but we can deduce it from the subject name
rc_data <- full_join(rc_data_temp, rc_input_csv, by=c("subject_id", "timepoint")) %>%
mutate(sex = as.factor(ifelse(grepl("F", subject_id), "female", "male")))
#
# Analyze each time point
#
# in this example those are the 9 time points, we'll show only for one time point
rc_data_p65 <- subset(rc_data, rc_data$timepoint == 65 & rc_data$fwhm == 0.2)
# voxel analysis, absolute Jacobians, fwhm 0.2
lm_abs_jac_p65 <- mincLm(log_det_absolute_second_level ~ sex, rc_data_p65, mask="mask.mnc")
qvals_lm_abs_jac_p65 <- mincFDR(lm_abs_jac_p65, mask="mask.mnc")
common_avg <- mincGetVolume("[pipe_name]_nlin_common/nlin-3.mnc")
# results at 10% FDR
mincPlotSliceSeries(mincArray(common_avg), mincArray(lm_abs_jac_p65, "tvalue-sexmale"), anatLow=500, anatHigh=2500, low=attr(qvals_lm_abs_jac_p65, "thresholds")[3,3], high=10, begin=40, end=-40, legend="T-statistics", mfrow=c(6,5), symmetric = TRUE)
# voxel analysis, relative Jacobians, fwhm 0.2
lm_rel_jac_p65 <- mincLm(log_det_relative_second_level ~ sex, rc_data_p65, mask="mask.mnc")
qvals_lm_rel_jac_p65 <- mincFDR(lm_rel_jac_p65, mask="mask.mnc")
# results at 15% FDR
mincPlotSliceSeries(mincArray(common_avg), mincArray(lm_rel_jac_p65, "tvalue-sexmale"), anatLow=500, anatHigh=2500, low=attr(qvals_lm_rel_jac_p65, "thresholds")[4,3], high=10, begin=40, end=-40, legend="T-statistics", mfrow=c(6,5), symmetric = TRUE)
# structurewise analysis (NB: use unblurred absolute Jacobians!)
rc_data_p65_no_blur <- subset(rc_data, rc_data$timepoint == 65 & rc_data$fwhm == 0.0)
# at present the registration chain code doesn't do any atlas generation; see "atlas generation" above
vols_p65 <- anatGetAll(rc_data_p65_no_blur$log_det_absolute_second_level, atlas = "your_atlas.mnc", defs = "your_atlas_mapping_of_labels.csv")
vols_p65_combined <- anatCombineStructures(vols_p65, defs="your_atlas_mapping_of_labels.csv")
anatlm_out <- anatLm(~ sex, rc_data_p65_no_blur, vols_p65_combined)
anatFDR(anatlm_out)
anatlm_out[order(abs(anatlm_out[,"tvalue-sexmale"])),]
Here we give an example for a single structure only...
rc_data_long_no_blur <- subset(rc_data, rc_data$fwhm == 0.0)
anat <- anatGetAll(rc_data_long_no_blur$log_det_absolute_second_level, atlas=your_atlas, method="jacobians", defs=your_label_defs)
anat_combined <- anatCombineStructures(anat, defs=your_label_defs)
rc_data_long_no_blur$bnst <- anat_combined[, "bed nucleus of stria terminalis"]
ggplot(data=rc_data_long_no_blur, aes(x=timepoint, y=bnst, group=subject_id, color=sex)) + geom_point(alpha=0.3)+geom_line(alpha=0.3)+stat_smooth(aes(group=sex))