Embryo scan preprocessing
This set of instructions is aimed at preprocessing embryo scans in minc format (acquired at MICe) in order to analyze differences in brain volume.
Step 1: Orient your scans
Flip the raw data in the yz direction to re-orient it in the same direction as the E18 atlas, which can be found here /opt/quarantine/resources/c57_E18_atlas/E18.0.mnc.
You can do so with the following command:
for file in ../raw_data/*.mnc; do volmash -swap yz $file $(basename $file .mnc)_swapped.mnc; done
Step 2: Manual corrections Some of your scans may violate assumptions made by the preprocessing algorithms, such as orientation or "extra" tissue/artefacts outside of the brain. Before we can automatically process, we may need to perform some manual corrections. To do so, we first create an Otsu mask that separates the embryo's body from the rest of the image with the following steps:
#In the directory with your swapped images
#1. get rid of negative intensities
for file in ./*.mnc; do mincmath -clobber -clamp -const2 0 1.79769e+308 $file ../cleaned/$(basename $file .mnc)_clamped.mnc;done
#2. Threshold the image to creat the otsu mask
module load ANTs
for file in /data/chamal/projects/lani/embryos/cleaned/*; do ThresholdImage 3 $file ./$(basename $file .mnc)_otsu.mnc Otsu 2;done
#3. Assign voxels within the mask 1 and outside the mask 0
for file in /data/chamal/projects/lani/embryos/cleaned/*otsu*; do ThresholdImage 3 $file ./$(basename $file .mnc)_otsu2.mnc 1 Inf 1 0; done
#4.Find the dimensions (?) of the largest component from the otsu
for file in ./*otsu2*; do ImageMath 3 ./$(basename $file otsu_otsu2.mnc)_otsu3.mnc GetLargestComponent $file; done
#5. Fill in any holes
for file in ./*otsu3*; do ImageMath 3 ./$(basename $file _otsu3.mnc)otsu4.mnc FillHoles $file 2; done
Next we will use Display to overlay the otsu mask (label) on top of the image to see if it masks the roi properly.
Display [swapped_image.mnc] -label [corresponding_otsu4.mnc]
Here is an example of an image with the extra tissue in the label:
Using the following process, we can fix the mask so that it looks more like this image:
Please note, this is a rough mask to help our preprocessing, it does not need to be perfect to the voxel level.
As usualy, the Display command will pop open the image as well as a on-screen keyboard. You can either use your physical keyboard or that virtual one.
First, hit space to make sure you pop to the top menu. Then hit D Colour Coding D Gray Scale, to turn the background image grey.
Hit space Pop menu to return to the main menu.
Hit F Segmenting. F XY Radius: 3 . Enter a value for the brush size I like to use 0.2 for fine segmenting.
Hit G Out-plane Rad: 0. Enter a value for the brush size in the out-of-plane dimension. I like to use 0.1.
To zoom in on a plane, scroll.
To move the center of the image, hold shift and left click and drag the image. Pretty intuitive!
To add more voxels to the label, right click your mouse.
To erase voxels from the label, hold shift and left click the mouse.
If you mess up, you can use ctrl - z.
Use the + and - keys to move through the image.
When you are ready to save, use space to pop to the main menu. Hit T File. Hit W Save Labels.mnc Choose a name/file location to save your label. Good job!
After you have generated all of your fixed labels (only took a couple of minutes, right? ;D) it's time to apply them to our original images. The label comprises a value of 1 for regions we want to include and a 0 for regions we don't want to include, so to apply this label to an image, we just need to multiply the original image and the label file since all voxels outside of the ROI will be set to 0. There are different tools we can use to do this. Below is one option to loop through all mncs in a folder named swapped and to multiply that file by the corresponding label file in a folder named labels before storing the output in a folder named cropped_to_label.
for file in ./swapped/*.mnc;
do mincmath -mult $file ./labels/$(basename $file _sept2014_dist_corr_swapped.mnc)*.mnc ./cropped_to_label/$(basename $file _sept2014_dist_corr_swapped.mnc)_labelled.mnc;
done
If we look at the output for the same file that we visualized earlier we see this:
Step 3: Preprocess the scans
Save the following scrip, written by Dr. Gabriel A Devenyi, as, for example, embryo-preprocessing.sh:
#!/bin/bash
set -euo pipefail
set -x
tmpdir=$(mktemp -d)
input=$1
output=$2
mincmath -clamp -const2 0 $(mincstats -max -quiet ${input}) ${input} ${tmpdir}/clamp.mnc
volcentre -com ${tmpdir}/clamp.mnc ${tmpdir}/centered.mnc
minc_anlm --rician --mt $(nproc) ${tmpdir}/centered.mnc ${tmpdir}/centered.denoise.mnc
mv -f ${tmpdir}/centered.denoise.mnc ${tmpdir}/centered.mnc
ThresholdImage 3 ${tmpdir}/centered.mnc ${tmpdir}/thresholdmask.mnc Otsu 1
ImageMath 3 ${tmpdir}/thresholdmask.mnc GetLargestComponent ${tmpdir}/thresholdmask.mnc
ImageMath 3 ${tmpdir}/thresholdmask.mnc FillHoles ${tmpdir}/thresholdmask.mnc 2
minccalc -unsigned -byte -expression '1' ${tmpdir}/centered.mnc ${tmpdir}/fullmask.mnc
N4BiasFieldCorrection -d 3 -s 8 -i ${tmpdir}/centered.mnc -b [20] -c [300x300x300x300,1e-6] \
-w ${tmpdir}/thresholdmask.mnc -x ${tmpdir}/fullmask.mnc -o ${tmpdir}/N4.mnc --verbose
ThresholdImage 3 ${tmpdir}/N4.mnc ${tmpdir}/thresholdmask.mnc Otsu 1
ImageMath 3 ${tmpdir}/thresholdmask.mnc GetLargestComponent ${tmpdir}/thresholdmask.mnc
ImageMath 3 ${tmpdir}/thresholdmask.mnc FillHoles ${tmpdir}/thresholdmask.mnc 2
N4BiasFieldCorrection -d 3 -s 4 -i ${tmpdir}/centered.mnc -b [20] -c [300x300x300x300,1e-6] \
-w ${tmpdir}/thresholdmask.mnc -x ${tmpdir}/fullmask.mnc -o ${tmpdir}/N4.mnc --verbose
ThresholdImage 3 ${tmpdir}/N4.mnc ${tmpdir}/thresholdmask.mnc Otsu 1
ImageMath 3 ${tmpdir}/thresholdmask.mnc GetLargestComponent ${tmpdir}/thresholdmask.mnc
ImageMath 3 ${tmpdir}/thresholdmask.mnc FillHoles ${tmpdir}/thresholdmask.mnc 2
ImageMath 3 ${tmpdir}/N4.mnc TruncateImageIntensity ${tmpdir}/N4.mnc 0.001 0.9995 2048 ${tmpdir}/thresholdmask.mnc
ImageMath 3 ${tmpdir}/N4.mnc RescaleImage ${tmpdir}/N4.mnc 0 65535
ImageMath 3 ${tmpdir}/N4.mnc PadImage ${tmpdir}/N4.mnc 50
antsApplyTransforms -d 3 -i ${tmpdir}/thresholdmask.mnc -o ${tmpdir}/thresholdmask.mnc -r ${tmpdir}/N4.mnc -n GenericLabel
iMath 3 ${tmpdir}/thresholdmask.mnc MD ${tmpdir}/thresholdmask.mnc 4 1 ball 1
ImageMath 3 ${tmpdir}/N4.mnc m ${tmpdir}/N4.mnc ${tmpdir}/thresholdmask.mnc
ExtractRegionFromImageByMask 3 ${tmpdir}/N4.mnc ${output} ${tmpdir}/thresholdmask.mnc 1 20
rm -rf ${tmpdir}
You can then run it with the following line:
for file in *.mnc; do echo ./embryo-preprocessing2.sh $file N4/$(basename $file swapped.mnc)preproc.mnc; done > joblist2
If using the command above, make a directory called N4 in your working directory to store the output folders.
To run it:
module load minc-toolkit minc-toolkit-extras qbatch
qbatch --ppj 8 joblist2
Step 4: Crop images to include just the brain
Now that you have your beautiful N4 corrected images, we will create an average of them so that we can crop out just the brain.
AverageImages 3 bulkaverage.mnc 1 *mnc
Next, you need to open this in register and create 8 tags to make a cube around the brain. Save these tags as a .tag file. Then you want to convert your .tag file to .mnc file with the following command (this will convert your .tag file to .mnc file, like the bulkaverage.mnc file):
tagtominc bulkaverage.mnc crop_box.tag tags.mnc
Now you can use this minc file to crop your average with the command below. This will apply the box that you created (tags.mnc) to crop the bulkaverage (input), and give you your croppedaverage.mnc (output) file.
autocrop -bbox tags.mnc bulkaverage.mnc croppedaverage.mnc
In my case, some of the scans were cut off with that box, so I ran the command below. the -extend option allows you to expand your bounding box to include ore of the scan. You can then run you joblist (as described above).
for file in *preproc.mnc; do echo autocrop -bbox tags.mnc $file ../cropped_inputs/$(basename $file .mnc)_cropped.mnc -extend 0v,20v 0v,0v 30v,30v -clobber; done > joblist_crop
**Optional final step:
In my case, some of my embryos had extra pieces of tissue outside the head that were causing some of the registrations to fail. For that, I used Display to create a paint a label that covered the extra tissue (for each individual scan that had problematic tissue; by going into the Segmenting panel (F) in Display you can reduce the size of your paintbrush by changing the dimensions of the XY radius (F), and then use your right click button of your mouse to shade in the area you want to remove. If you make a mistake you can undo or click shift plus right click and erase)
To invert the labels so that you can delete the extra bits outside the embryo:
for file in *label.mnc; do ImageMath 3 $(basename $file .mnc)_invert.mnc Neg $file; done
To multiple the embryo scans by the inverted masks to delete the extra bits:
for file in path_to_labels_for_cleanup/*invert.mnc; do ImageMath 3 directory_for_cleaned_up_images/$(basename $file _label_invert.mnc)_cleaned.mnc m $file $(basename $file _label_invert.mnc).mnc; done