Introduction
The twolevel_model_building.py code was written for the following use case: all subjects in your data set are scanned multiple times, but in contrast to the type of longitudinal data used for the longitudinal registration tools (registration_chain.py), all timepoints for a given subject can be registered together. This is done using iterative group-wise registration to create a subject-specific average. All of these averages are then registered together, again using a iterative group-wise procedure, to create a population average.
Typical run
The input files to the pipeline are read from a comma separated file containing the (string) columns 'file' and 'group' (and possibly others); for instance:
> cat TwoLevelFiles.csv file,group subject_1_time_point_1.mnc,1 subject_1_time_point_2.mnc,1 subject_1_time_point_3.mnc,1 subject_1_time_point_4.mnc,1 subject_1_time_point_5.mnc,1 subject_2_time_point_1.mnc,2 subject_2_time_point_2.mnc,2 subject_2_time_point_3.mnc,2 subject_3_time_point_1.mnc,3 subject_3_time_point_2.mnc,3 subject_3_time_point_3.mnc,3 subject_1_time_point_4.mnc,3 subject_4_time_point_1.mnc,4 subject_4_time_point_2.mnc,4 subject_4_time_point_3.mnc,4
A sample call of the program could be:
twolevel_model_building.py \ --num-executors=600 \ --init-model=/hpf/largeprojects/MICe/tools/initial-models/Pydpiper-MEMRI-90micron-saddle-july-2015/p65_MEMRI_mouse_brain.mnc \ --pipeline-name=test \ --lsq12-protocol=/hpf/largeprojects/MICe/tools/protocols/linear/Pydpiper_testing_default_lsq12.csv \ --nlin-protocol=/hpf/largeprojects/MICe/tools/protocols/nonlinear/Pydpiper_mincANTS_standard_90_micron.pl \ --registration-method=ANTS \ --output-dir=test \ --no-run-maget \ --default-job-mem=8 \ --maget-no-mask \ --csv-file=TwoLevelFiles.csv
What the output looks like
The first level registrations will be stored in the directory:
twolevel_brains_firstlevel # which will contain directories belonging to the iterative build models for each of the input subjects: subject_1_time_point_1_lsq12 subject_1_time_point_1_lsq6 subject_1_time_point_1_nlin subject_1_time_point_1_processed subject_2_time_point_1_lsq12 subject_2_time_point_1_lsq6 ...
After those have been created, the second level registration (iterative build model on all the final non linear averages from the first level) will be run:
twolevel_brains_secondlevel # which will contain: second_level_lsq12 second_level_lsq6 second_level_nlin second_level_processed
Manual verification of the pipeline
Generating the following images will give you a sense of how well the registration worked.
#
# the first level registrations:
#
# set scale to 1 for human data, and to 20 for mouse data
scale=20;
for nlindir in *_first_level/*_processed;
do topleveldir=`dirname $nlindir`;
subjectbase=`basename $nlindir _nlin`;
for subj in ${topleveldir}/${subjectbase}/*/resampled/*N_I_lsq6_lsq12_and_nlin-resampled.mnc;
do subjbase=`basename $subj .mnc`;
mincpik -clobber -scale $scale --triplanar $subj /tmp/${subjectbase}_sample_${subjbase}.png;
done;
mincpik -clobber -scale $scale --triplanar ${nlindir}/*/*nlin-3.mnc /tmp/${subjectbase}_avg.png;
montage -geometry +2+2 /tmp/${subjectbase}_avg.png /tmp/${subjectbase}_sample_*png ${topleveldir}/${subjectbase}_QC.png;
done
# view the image from the command line:
eog *_first_level/*_QC.png
#
# the second level registrations:
#
for secondlevelnlin in *_second_level/*_nlin;
do secondleveltopdir=`dirname $secondlevelnlin`;
for subjavg in ${secondleveltopdir}/second_level_processed/*/resampled/*final-nlin.mnc;
do subjavgbase=`basename $subjavg .mnc`;
mincpik -clobber -scale $scale --triplanar $subjavg /tmp/second_level_sample_${subjavgbase}.png;
done
mincpik -clobber -scale $scale --triplanar *_secondlevel/*_nlin/*nlin-3.mnc /tmp/second_level_avg.png;
montage -geometry +2+2 /tmp/second_level_avg.png /tmp/second_level_sample_*.png ${secondleveltopdir}/secondlevel_QC.png
done
eog *_secondlevel/*_QC.png
Stats volumes
The original intended purpose of the twolevel_model_building code can be explained by the following example: say you have a drug and you want to determine the effects the drug has on brain growth. The brain growth can be determined from the deformations (Jacobian determinant of the deformations) from the first level registrations. In order to compare brain growth between your subjects, these Jacobian determinants need to live in the same space. This is what the second level registrations are for. The Jacobian determinants from the first level are resampled into the second level final average space, and this is what you will do you statistics on:
# All stats files reside in the _firstlevel directories. For each subject you will have an _processed directory, with directories for each time point. The stats files you will find have the following naming: # absolute Jacobians from first level ...-final-nlin_with_additional_inverted_absolute_log_determinant_fwhm0.1.mnc ...-final-nlin_with_additional_inverted_absolute_log_determinant_fwhm0.2.mnc ... # relative Jacobians from first level (overall scaling has been taken out) ...-final-nlin_with_additional_inverted_pure_nlin_relative_log_determinant_fwhm0.1.mnc ...-final-nlin_with_additional_inverted_pure_nlin_relative_log_determinant_fwhm0.2.mnc and # absolute Jacobians from first level resampled into the final nlin space of the second level ...-final-nlin_with_additional_inverted_absolute_log_determinant_fwhm0.1_common.mnc ...-final-nlin_with_additional_inverted_absolute_log_determinant_fwhm0.2_common.mnc ... # relative Jacobians from first level resampled into the final nlin space of the second level ...-final-nlin_with_additional_inverted_pure_nlin_relative_log_determinant_fwhm0.1_common.mnc ...-final-nlin_with_additional_inverted_pure_nlin_relative_log_determinant_fwhm0.2_common.mnc
There is another way to link all the input files together. This is by concatenating the transformations from a subject to the first level average, and concatenating that with the transformation from that first level average to the second level average and determining the Jacobian determinants based on those transformations. This is currently not being done by the twolevel_model_building.py pipeline.
Generating stats originating in the second level average
If you are using pydpiper version 1.17 or earlier you can generate the statistic files originating in the second level average as follows. The code below has some comments to indicate what's happening. In the second code block all the commands are submitted to sge, to have the processing be done on the farms.
# Assumption: the pipeline basename in this example is: twolevel_pipe
# Files that will be generated for each input file:
#
# Absolute log determinants from the second level average:
# */stats-volumes/*-transformation-from-lsq6-to-second-level-final-average_inverse_grid_log_determinant.mnc
# */stats-volumes/*-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_log_determinant.mnc
#
# Relative log determinants from the second level average:
# */stats-volumes/*-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_log_determinant.mnc
# */stats-volumes/*-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_log_determinant.mnc
# 1) concatenate the full transformation from the first level with the full transformation from the second level:
pipeline_base=twolevel_pipe
for file in ${pipeline_base}_firstlevel/*_processed;
do first_level_group=`basename $file _processed`;
for subject_sub_dir in $file/*;
do subject=`basename $subject_sub_dir`;
# concatenate transforms:
xfmconcat ${subject_sub_dir}/transforms/${subject}-final-nlin_with_additional.xfm ${pipeline_base}_secondlevel/second_level_processed/${first_level_group}-nlin-3/transforms/${first_level_group}-nlin-3-final-nlin_with_additional.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average.xfm
xfminvert ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse.xfm
minc_displacement ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse.xfm ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid.mnc
#
# absolute Jacobians
#
mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant_temp.mnc
mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant.mnc
mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_log_determinant.mnc
#
# absolute Jacobians with 0.2mm Gaussian blur (change if you want)
#
smooth_vector --filter --fwhm=0.2 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2.mnc
mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant_temp.mnc
mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant.mnc
mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_log_determinant.mnc
#
# relative Jacobians ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc
#
lin_from_nlin ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_linear_part.xfm
xfmconcat ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_linear_part.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only.xfm
minc_displacement ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only.xfm ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid.mnc
mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant_temp.mnc
mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant.mnc
mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_log_determinant.mnc
#
# relative Jacobians with 0.2mm Gaussian blur
#
smooth_vector --filter --fwhm=0.2 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2.mnc
mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant_temp.mnc
mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant.mnc
mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_log_determinant.mnc
done
done
Submit the whole things to the farms:
# Assumption: the pipeline basename in this example is: twolevel_pipe
# 1) concatenate the full transformation from the first level with the full transformation from the second level:
pipeline_base=twolevel_pipe
for file in ${pipeline_base}_firstlevel/*_processed;
do first_level_group=`basename $file _processed`;
for subject_sub_dir in $file/*;
do subject=`basename $subject_sub_dir`;
sge_batch -l vf=6G xfmconcat ${subject_sub_dir}/transforms/${subject}-final-nlin_with_additional.xfm ${pipeline_base}_secondlevel/second_level_processed/${first_level_group}-nlin-3/transforms/${first_level_group}-nlin-3-final-nlin_with_additional.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average.xfm \; xfminvert ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse.xfm \; minc_displacement ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse.xfm ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid.mnc \; mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant_temp.mnc \; mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant.mnc \; mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_log_determinant.mnc \; smooth_vector --filter --fwhm=0.2 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2.mnc \; mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant_temp.mnc \; mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant.mnc \; mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_grid_fwhm0.2_log_determinant.mnc \; lin_from_nlin ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_linear_part.xfm \; xfmconcat ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_linear_part.xfm ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only.xfm \; minc_displacement ${pipeline_base}_secondlevel/second_level_nlin/second_level-nlin-3.mnc ${subject_sub_dir}/transforms/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only.xfm ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid.mnc \; mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant_temp.mnc \; mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant.mnc \; mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_log_determinant.mnc \; smooth_vector --filter --fwhm=0.2 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2.mnc \; mincblob -det ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant_temp.mnc \; mincmath -add -const 1 ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant_temp.mnc ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant.mnc \; mincmath -log ${subject_sub_dir}/tmp/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_determinant.mnc ${subject_sub_dir}/stats-volumes/${subject}-transformation-from-lsq6-to-second-level-final-average_inverse_nonlinear_only_grid_fwhm_0.2_log_determinant.mnc
done
done