Child pages
  • Two level registration (twolevel_model_building)
Skip to end of metadata
Go to start of metadata

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