A landmark-free morphometrics pipeline for high-resolution phenotyping: application to a mouse model of Down syndrome

ABSTRACT Characterising phenotypes often requires quantification of anatomical shape. Quantitative shape comparison (morphometrics) traditionally uses manually located landmarks and is limited by landmark number and operator accuracy. Here, we apply a landmark-free method to characterise the craniofacial skeletal phenotype of the Dp1Tyb mouse model of Down syndrome and a population of the Diversity Outbred (DO) mouse model, comparing it with a landmark-based approach. We identified cranial dysmorphologies in Dp1Tyb mice, especially smaller size and brachycephaly (front-back shortening), homologous to the human phenotype. Shape variation in the DO mice was partly attributable to allometry (size-dependent shape variation) and sexual dimorphism. The landmark-free method performed as well as, or better than, the landmark-based method but was less labour-intensive, required less user training and, uniquely, enabled fine mapping of local differences as planar expansion or shrinkage. Its higher resolution pinpointed reductions in interior mid-snout structures and occipital bones in both the models that were not otherwise apparent. We propose that this landmark-free pipeline could make morphometrics widely accessible beyond its traditional niches in zoology and palaeontology, especially in characterising developmental mutant phenotypes.


Cranial Interior morph stretch heatmap
Cranial morph, showing the interior of the skull, between the mean WT and Dp1Tyb shapes, where size has been regressed out. Red and Blue colour change indicate expansion and contraction respectively as a percentage. Derived from landmark-free analysis.

Cranial mesh points stretch heatmap
Landmark-free-derived heatmap of cranium local stretch with no scaling. The mesh is represented here as all the points of the mesh rather than a rendered surface. Red and Blue colour change indicate expansion and contraction respectively. Downloading and viewing in a paused video by moving the slider is recommended. Development: doi:10.1242/dev.188631: Supplementary information Movie 9.

Mandible morph stretch heatmap
Mandible morph between the mean WT and Dp1Tyb shapes, where size has been regressed out. Red and Blue colour change indicate expansion and contraction respectively as a percentage. Derived from landmark-free analysis.

Cranial morph stretch heatmap
Cranial morph between the mean WT and Dp1Tyb shapes, where size has been regressed out. Red and Blue colour change indicate expansion and contraction respectively as a percentage. Derived from landmark-free analysis. Development: doi:10.1242/dev.188631  Planar view of mouse head captured by µCT. A Image pre-processing. B Image postprocessing. Yellow arrow indicates cartilaginous element (nasal turbinate). Development: doi:10.1242/dev.188631: Supplementary information : doi:10.1242/dev.188631: Supplementary information Figure S5. Validation of landmark-free morphometry in a sample of Diversity Outcross Mice by quantification of sexual dimorphism. A. Canonical Variance Analysis by sex; i) distributions of CV scores by method (unlike the DO group, these specimens are wellseparated into two distinct groups); ii) Bland-Altman plot for CV scores for sexual dimorphism. Here, the deviations of each method from the manual landmarking-based estimated is plotted against the manual landmarking estimate; iii) variance components for sexual dimorphism and correlations among the CV scores across methods. B. 3D Morphs and heatmaps showing the shape variation associated with sexual dimorphism as estimated by each method. i) GM = manual landmark-based geometric morphometics; ii) GP OPT = automated landmark-based geometric morphometics; iii) LF 2808 and iv) LF 560 = landmark-free method with 2808 or 560 control points respectively. "Min x2" and "Max x2" label the deformations corresponding to the lower and higher ends of the heat map colour scales, exaggerated two-fold to more easily visualise the shape differences.  Table S1.

Development • Supplementary information
Comparison of normalised centroid sizes of crania and mandibles from WT and Dp1Tyb mice determined using landmark-based or landmark-free methods.

Landmark-free Morphometrics Pipeline Description
This Appendix is intended as an overview of the landmark-free pipeline for a relative nonspecialist. Detailed dcoumentation is embedded with the pipeline code published on GitLab at https://gitlab.com/ntoussaint/landmark-free-morphometry. This Appendix is organised into five sections as follows: A. Pipeline Dependencies B. µCT Pre-processing C. Mesh Alignment D. Atlas Construction E. Shape Statistics

A. Pipeline dependencies and notes on initialisation
All packages and software listed below are required. The pipeline is written in Python and C++ and requires a suitable environment to run both. We used the browser-based Jupyter environment (https://jupyter.org/). Specific requirements for each segment of the pipeline are listed at the appropriate page within the pipeline code but the complete list is as follows: The pipeline code is downloaded as landmark-free-morphometry-master.zip or ~.tar. After decompressing the package and moving to its directory within the Terminal app (Apple, Linux or equivalent in Windows). Jupyter is started using the command ipython notebook which opens the pipeline in the default system Web browser. The pipeline is presented as a series of virtual "pages" which contain "cells" (sections) of two types. One type defines a function, the other applies those functions to the data. Cells can be run individually or as an automatic sequence (see Jupyter documentation for details). An example of a function-defining cell is below:

B. µCT Pre-processing
This section describes the code used to extract and produce a mesh from an initial µCT image.
The following cells can be run with default values to define the required functions: Binary Segmentation -Defines functions used in the object extraction Parcellation -Defines functions required for parcellation of binary image Mesh Extraction -Defines function used in generating a mesh from a binary image Mirror mesh file -Defines function used to mirror the mesh file Display help functions -Defines functions used to display images within the notebook Generic Imports -Imports Numpy Using ITK-Snap save µCT images in the NIFTI (.nii.gz) format.
The following five cells need to be edited to change parameters and filenames before running (N.B. annotations describe the sections that can or need to be changed, along with a brief description of inputs, outputs and function/s of the cell).

Load data
Loads in the image file to be processed.

Object extraction
The cell uses the functions defined in object_extraction to extract the object of interest from the input image and apply morphological functions to remove any noise and close unwanted holes.
Defines current working directory. Set the folder in which outputs are stored. This should not need to be changed as the 'data' and 'preprocessing' folders are already within the file along with the scripts.
In this example uCT-1.nii.gz is the initial µCT image to be processed. This should be replaced with the filename desired.
This tells the script where to look for the file to be processed. In this case within the data folder.
Check output file using ITK-Snap to validate successful extraction of the region of interest.
If the regions has been appropriately binarized it should be saved as .mha file.
Name of output file, a binarized image.

Adjustable parameters
Threshold -an appropriate value should be selected to extract the region of interest from the original µCT image. Closing_kernel and opening_kernel -parameters used to fill holes in the binary image. If none is required set to 0 Development: doi:10.1242/dev.188631: Supplementary information Parcellation This cell uses the functions defined in parcellation to 'parcellate' the image into various regions according to the density of the bone. This allows the separation of mandible (densest bone) from the rest of the skull. If no region separation is required skip this step and move to the Mesh extraction step: The output is a binary image with regions ordered by bone density. Use ITK-Snap to visualise these labels. The example on the right shows these regions, here label one is the red region. Note the label/s for region of interest (ROI).

Object selection
Input the labels of relevance to separate the ROI from the rest of the subject.
Check the output in ITK-snap opening it as a segmentation file. Any noise or irrelevant regions can be removed using the paint brush tool.
Define the name for file to be segmented i.e. the file from Object extraction.
Watershed level should be .22 for full resolution image. ( .1 is for the example data provided).
Define the file name for the selected region of interest. This is the output file from region parcellation.
Input labels of regions to be extracted from the region parcellation step prior.
If only the first label is required, input 1 for both the lower threshold (l_threshold) and upper threshold (u_threshold). If the first 2 labels are desired set the values to 1 and 2 respectively. If the whole image is desired input the lowest and highest label numbers.

Extract mesh
This step uses functions defined in Mesh Extraction to extract a surface mesh from the region of interest selected prior. The meshes generated are in the .vtk format.
Mesh files produced are large in size and increase computational time of the Atlasing step to unrealistic durations. To overcome this, meshes were decimated using meshlab.
Save meshes in the .stl (stereolithography) format using paraview. Open the mesh in meshlab and apply quadric edge collapse decimation fi and enter the percentage reduction desired (authors used 0.15).
Once meshes for all subjects are generated move to the next step, Mesh Alignment.
Use paraview again to save the newly decimated mesh as a .vtk file and ensure it is saved in ASCII format Output mesh file name.

C. Mesh Alignment
This section describes the process of aligning the meshes produced in the previous section. Initial coarse alignment is required to allow the Atlasing process (below) to work. This first requires manually placing a small number of landmarks using Mitk-Workbench, minimally three (e.g. for a hemi-mandible) or more usually three symmetrically placed pairs. Landmark file names should be identical to the mesh file name e.g. Mandible_1.vtk, with its landmark file, Mandible_1.mps. Below is shows landmarks used to align the mandible.
Once complete move landmark files (.mps) and mesh files (.vtk) into Alignment folder.
The following are run with default values:

Read landmarks -Defines functions used to import landmark files Numpy -vtk help functions
-Define functions used convert landmark files so they can be used to align meshes -Functions used to extract centroid and centroid size Procrustes alignment -Defines functions used in alignment o Rigid -rigid body alignment o Similarity -rigid body plus scaling alignment Load data Loads meshes and associated landmark files ready for alignment, Development: doi:10.1242/dev.188631: Supplementary information Mesh Alignment Uses functions defined in Procrustes alignment to align meshes. Two modes can be selected, Similarity or Rigid body. "Similarity" aligns meshes whilst regressing out the size. "Rigid body" aligns meshes but maintains size differences: Mesh alignment fidelity should be checked in Paraview.

D. Atlas Construction
This section describes the application of the shape averaging and comparison (atlasing) processes. Within the atlasing folder (included in the pipeline package) there are three files adjustable files: Model.xml describes parameters for deformation: -kernel-width: Order of magnitude of displacements for the template -kernel-width: Order of magnitude for the subjects' displacements Here, define the folder in which the pipeline will search for inputs and store outputs. By default, it is set to the alignment folder within data.
Choose mode here to 'Similarity' or 'Rigid body' Adds a suffix to aligned mesh files so they can easily be distinguished from the original unaligned files. Development: doi:10.1242/dev.188631: Supplementary information Optimization_parameters.xml describes parameters that drive the minimization procedure: -max-iterations: Stop criterion for iterations, usually ~150 -number-of-threads: Number of threads to use, usually equal to the number of subjects (should not exceed the number of threads available on the system) Data_set.xml defines the files to be atlased and should be structured as follows: From the alignment folder transfer aligned meshes (files with suffix _r) and initialtemplate.vtk into the atlas folder.

Deformetrica
The first cell loads the Deformetrica software. If the software is not found an error is returned.

Launch Simulation
Launches the atlasing software using the files and parameters defined by model.xml, optimisation_parameters.xml and data_set.xml. Deformetrica.log, within the output folder, contains information on the progress of the atlas, including the number of control points and the current iteration.

E. Shape statistics
This section describes the statistical analysis of the outputs from the atlas. Atlasing generates several output files of which the following three are used for statistical analysis: -Atlas_controlpoints.txt o Control points of the atlas -Atlas_momenta.txt o Momentum vectors of the atlas -Atlas_initial_template.txt o Average mesh of the population These files are must be moved into the shape statistics folder. In the same folder, a userpopulated file data.csv defines the names and types of the subjects. data.csv example file In this example in GroupId, 1 refers to a mutant and -1 to WT, and in Gender, 1 Male and -1 female.
The following cells are the run with default values: Imports -Loads relevant packages required for this section of the pipeline Deformetrica -Loads deformetrica ready to be used (required for creation of morphs) Load data -Loads output files stated previously to be analysed Define population's groups -Loads data.csv Development: doi:10.1242/dev.188631: Supplementary information Subgroup definitions -Assigns names to groups of specimens defined by one or more parameters in data.csv.
[Default groups are: WTf (Wildtype female, defined by: group ID = -1, gender = -1, CRf is defined by: group id = 1, gender =-1), WTm (Wildtype Male), CRf (Carrier female) and CRm (Carrier Male)] However, the user should change the names and definitions of subgroups to suit their needs.

Kernel Principal Component Anaylsis (kPCA)
This cell runs a principal component analysis (PCA) and produces a classification score (as a percent) and a p value for said score. As well as defining group IDs (by default group one > 0 and group two < 0 using values set in data.csv) and calculating the mean Principal Component score. Development: doi:10.1242/dev.188631: Supplementary information Eigenvalues/Eigenvectors Produces Eigenvalues and eigen vectors. Providing data on how much variability each principle component describes.

Write PCA points on disk
Produces an excel file, kpca.csv containing the PC values which can exported to produce PCA graphs.

Momenta projection
Produces momenta that describes project between the mean of the whole population to the means of the previously defined populations.

Shoot mean shape between groups
Uses the momenta projection to 'shoot' the shape towards the means of the two groups.

Shooting outputs
Saves the shot meshes produced in the previous step in two locations: data/shapestatistics/shooting/forward data/shapestatistics/shooting/backward Cross validation score and standard deviation Cross validation score and p value of permutation test Development: doi:10.1242/dev.188631: Supplementary information They correspond to the mean shape (Atlas_initial_template.vtk) displaced towards the mean of group one (backward) and the mean of group two (forward). Each file now contains several meshes. The first mesh in the folder is the average shape of the whole population and the last mesh the mean shape of the subgroup. Each mesh in-between represents a step of deformation moving from one to the other.

Animation of deformations between subgroups
Reorders the mesh files produced in Shooting outputs to generate an animation of cyclic shape change starting with the mean of group one to the mean of group two and back again. The files for the cyclic deformation animation are saved in: data/shapestatistics/shooting/combined. Files should be opened in Paraview to visualise the meshes.
The deformation is viewed as a cycling animation by pressing the repeat and play buttons. Development: doi:10.1242/dev.188631: Supplementary information "Stretch" map (Local volume change map) between groups This step integrates data on displacement and stretch into the cyclic deformation produced in the previous step. Removing the # from in front of compute_point_displacements calculates displacement.
Removing the # from in front of compute_cell_surfaces_difference, calculates stretch.
Heat maps are viewed by importing the mesh files from the combined folder into Paraview. Select displacements or Absolute Volume from the colouring drop down menu and the morph was played as before to visualise the heat maps. Colour scale (look-up table) may need to be adjusted to visualise differences properly. This can be found by clicking the edit button found underneath the dropdown menu.

Generating videos
To save these morphs as videos (either with heat map colours or not). Go to file and save animation.
This will open up a file save window. Here the format of the video (.avi as default) and the name of the video can be chosen. This will bring up a second window in which a number of parameters for the video such as resolution and frames per second maybe adjusted.
Once happy with the parameters hit okay and the video will be produced. The video will capture everything in the rendering window including colours scale bars. To edit the appearance of the scale bar, click the edit colour bar icon.
This will bring up a menu in which various settings of the colour scale bar may be changed.
To remove the scale bar entirely hit the colour bar icon. Development: doi:10.1242/dev.188631: Supplementary information