[Free|Pet]Surfer by Example

Init FreeSurfer env

At the NIH

• Connect to Biowulf via NoMachine
• Launch a terminal on Biowulf
• Run:

sinteractive --cpus-per-task 8 --mem 40G
module load freesurfer

In General

Point the environment variable $FREESURFER_HOME to the installation directory

export FREESURFER_HOME=/home/paul/lcn/freesurfer-7.3.2

Init the FreeSurfer environment

source $FREESURFER_HOME/SetUpFreeSurfer.sh

Add the following to ~.bashrc (Optional)

export FREESURFER_HOME=/home/paul/lcn/freesurfer-7.3.2
alias fs='source $FREESURFER_HOME/SetUpFreeSurfer.sh'

Set the Subject's directory

Create a directory to store FreeSurfer outputs

FreeSurfer calls this the "Subejct's Directory" and uses an environment variable $SUBJECTS_DIR to find it.

On Biowulf

export SUBJECTS_DIR=/data/$USER/fs-subjects
mkdir -p $SUBJECTS_DIR

On my machine

export SUBJECTS_DIR=/home/paul/lcn/20230519-freesurfer-training-nih/fs-subjects
mkdir -p $SUBJECTS_DIR

Define the location of the input data

We will use the openneuro dataset ds004230 as example data.

Let's use another environment variable to track the location of this dataset to help manage the differences between my machine and Biowulf

On Biowulf

export INPUT_DATADIR=/data/classes/freesurfercourse2023/pw

On my machine:

export INPUT_DATADIR=/home/paul/lcn/20230519-freesurfer-training-nih

Now that the following environment variables have been defined

• $FREESURFER_HOME
• $SUBJECTS_DIR
• $INPUT_DATADIR

We should be able to write all further commands only once, and have them work both on my machine and Biowulf.  We can use echo to verify these commands have been set

echo $FREESURFER_HOME
echo $SUBJECTS_DIR
echo $INPUT_DATADIR

Running recon-all

recon-all is the script that launches the main FreeSurfer processing pipeline.

A configuration file, global-expert-options.txt, if placed in $SUBJECTS_DIR, will override the default parameters of recon-all.

The recommended global-expert-options.txt files are located at:

• For FreeSurfer v7.3.2: $INPUT_DATADIR/global-expert-options-7.3.2.txt
• For FreeSurfer v7.4:  $INPUT_DATADIR/global-expert-options-7.4.txt

Biowulf is currently running FreeSurfer v7.3.2, so let's put that version of global-expert-options.txt in our $SUBJECTS_DIR

cp $INPUT_DATADIR/global-expert-options-7.3.2.txt $SUBJECTS_DIR/global-expert-options.txt

We are now ready to run recon-all, using the following flags:

• -i: location of input file
• -s: FreeSurfer subject name.  Output data will be in a folder under $SUBJECTS_DIR with this name
• -all: Run the whole pipeline

recon-all \
  -all \
  -i ${INPUT_DATADIR}/sub-PS19/ses-baselinebrain/anat/sub-PS19_ses-baselinebrain_T1w.nii.gz \
  -s sub-PS19-ses-baselinebrain

This command will take approx 3.5 hours and if everything worked the final few lines of output should look like:

Started at Sat May 20 10:25:14 EDT 2023 
Ended   at Sat May 20 13:56:59 EDT 2023
#@#%# recon-all-run-time-hours 3.529
recon-all -s sub-PS19-ses-baselinebrain finished without error at Sat May 20 13:56:59 EDT 2023
done

Since this command takes a while to run, let's just copy the output to our $SUBJECTS_DIR and pretend we've run it

cp -R ${INPUT_DATADIR}/fs-subjects/sub-PS19-ses-baselinebrain ${SUBJECTS_DIR}

More details on what recon-all is doing under the hood can be found here

Examining Outputs

Lets looks at the major directories

ls -l ${SUBJECTS_DIR}/sub-PS19-ses-baselinebrain

• label: Stores cortical labels
• mri: Stores MR volumetric data
• pet: Ignore for now, not generated by recon-all
• scripts: Details on how the data was generated
• stats: Statistics computed over the results (e.g cat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/stats/aseg.stats)
• surf: Stores Cortical surface meshes and surface overlays (like thickness, curvature, etc)
• tmp: Stores intermediates files (should be empty if the recon finished without errors)

Let's look at some outputs in FreeView

cd ${SUBJECTS_DIR}/sub-PS19-ses-baselinebrain
freeview \
  --volume ./mri/orig.mgz \
  --volume ./mri/aseg.mgz:opacity=0.5 \
  --surface ./surf/lh.white:edgecolor=yellow \
  --surface ./surf/rh.white:edgecolor=yellow \
  --surface ./surf/lh.pial:edgecolor=red:overlay=./surf/lh.thickness \
  --surface ./surf/rh.pial:edgecolor=red:annot=./label/rh.aparc.annot

Other useful commands to interrogate outputs:

mri_info ${SUBJECTS_DIR}/sub-PS19-ses-baselinebrain/mri/aseg.mgz
mris_info ${SUBJECTS_DIR}/sub-PS19-ses-baselinebrain/surf/lh.white

Convert from .mgz to .nii.gz

mri_convert \
  ${SUBJECTS_DIR}/sub-PS19-ses-baselinebrain/mri/aseg.mgz \
  ${SUBJECTS_DIR}/sub-PS19-ses-baselinebrain/mri/aseg.nii.gz

Petsurfer Step-by-Step

Reference

Segmentation for PVC

The conventional FreeSurfer segmentation (mri/aseg.mgz) doesn't label voxels outside of the brain (e.g. skull), so let's generate another more appropriate segmentation (you don't need to run this, it has already been run in the example data)

Running

gtmseg --s sub-PS19-ses-baselinebrain

Should create the file mri/gtmseg.mgz and produce the following output to the terminal if successful

...
Started at Sat 20 May 2023 04:03:45 PM EDT 
Ended   at Sat 20 May 2023 04:28:49 PM EDT
Gtmseg-Run-Time-Sec 1504
Gtmseg-Run-Time-Hours 0.42

Let's compare  mri/gtmseg.mgz with mri/aseg.mgz:

cd $SUBJECTS_DIR/sub-PS19-ses-baselinebrain
freeview \
  --volume mri/nu.mgz \
  --volume mri/aseg.mgz \
  --volume mri/gtmseg.mgz:colormap=lut

Examine PET data

freeview \
  --volume ${INPUT_DATADIR}/ds004230/sub-PS19/ses-baselinebrain/pet/sub-PS19_ses-baselinebrain_rec-DynTOF_pet.nii.gz

Create a pet folder for the subject

cd $SUBJECTS_DIR/sub-PS19-ses-baselinebrain
mkdir -p pet

Motion Correct PET data

mc-afni2 \
  --t $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz \
  --i $INPUT_DATADIR/ds004230/sub-PS19/ses-baselinebrain/pet/sub-PS19_ses-baselinebrain_rec-DynTOF_pet.nii.gz \
  --o $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/pet.mc.nii.gz \
  --mcdat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/motion-parameters.dat

Create PET template

By choosing an arbitrary frame:

cd $SUBJECTS_DIR/sub-PS19-ses-baselinebrain
mri_convert \
  ${INPUT_DATADIR}/ds004230/sub-PS19/ses-baselinebrain/pet/sub-PS19_ses-baselinebrain_rec-DynTOF_pet.nii.gz \
  --frame 5 \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz

output:

keeping frame 5
writing to ./pet/template.nii.gz...

By taking the mean across frames.  In this case it's better to use use motion corrected data as input:

cd $SUBJECTS_DIR/sub-PS19-ses-baselinebrain
mri_concat \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/pet.mc.nii.gz \
  --mean \
  --o $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz

View PET template and MR data

freeview \
  --hide-3d-slices \
  --volume $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/mri/orig.mgz \
  --volume $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz \
  --surface $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/surf/lh.white:edgecolor=yellow \
  --surface $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/surf/rh.white:edgecolor=yellow

Note that the data isn't aligned

Register PET image with the anatomical

If there is little uptake outside of the brain (eg, FDG and many receptor ligands, lets us this one):

mri_coreg \
  --s sub-PS19-ses-baselinebrain \
  --mov $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz \
  --reg $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.reg.lta

Viz the results of motion correction (before vs after motion correcting):

freeview \
  --volume $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/pet.mc.nii.gz \
  --volume $INPUT_DATADIR/ds004230/sub-PS19/ses-baselinebrain/pet/sub-PS19_ses-baselinebrain_rec-DynTOF_pet.nii.gz

verify:

freeview \
  --hide-3d-slices \
  --volume $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/mri/orig.mgz \
  --volume $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz:reg=$SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.reg.lta:colormap=grayscale \
  --surface $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/surf/lh.white:edgecolor=yellow \
  --surface $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/surf/rh.white:edgecolor=yellow

If there is lots of uptake outside the brain

mri_coreg \
  --s sub-PS19-ses-baselinebrain \
  --mov $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.nii.gz \
  --reg $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.reg.lta \
  --ref $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/mri/nu.mgz \
  --no-ref-mask

If you don't see anything, hold down shift and hold down the middle mouse button and drag a box around the red crosshair in one of the viewports.  This will automatically adjust the contrast

Partial Volume Correction and TAC generation

mri_gtmpvc \
  --i $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/pet.mc.nii.gz \
  --reg $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/template.reg.lta \
  --psf 6 \
  --seg $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/mri/gtmseg.mgz \
  --default-seg-merge \
  --auto-mask 1 .01 \
  --mgx .01 \
  --km-ref 8 47 \
  --km-hb 11 12 13 50 51 52 \
  --no-rescale \
  --o $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output

Where:

• --i: input data (motion corrected PET data)
• --o: the output folder
• --reg: registration from pet template to anatomical
• --psf: The full-width half max point spread function of your PET scanner, in mm (ask your PET physicist)
• --seg: The segmentation to use
• --default-seg-merge: Merge some smaller labels in seg together (eg, all the ventricular CSF segments will be merged into one ROI) before computing PVC
• --auto-mask FWHM threshold: create a mask so that irrelevant voxels outside of the are not analyzed
•  --mgx GM-threshold: Run "extended" Muller-Gartner analysis for all voxels that contain at least 1% grey matter
• --km-ref 8 47: Specifies the label values in seg to treat as the reference region
• --km-hb 11 12 13 50 51 52: Specifies the label values in seg to treat as high binding region

Values for --km-ref/--km-hb can be found via:

• cat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/mri/gtmseg.ctab
• FreeSurferColorLUT.txt
• freeview $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/mri/gtmseg.mgz:colormap=lut

Examine outputs:

cd $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output
ls -l
find . -maxdepth 1 -name '*.nii.gz' -exec sh -c 'echo {}; mri_info --dim {}; echo' \;

Some files (e.g. gtm.nii.gz) will have vox dims of 100 1 1 33.  These files are not viewable by FreeView or other viewers, but are stored as nifti to make them easy to use with tools like mri_glmfit

Outputs:

• km.ref.tac.dat: Reference TAC.  Average TAC computed across all voxels with labels specified by --km-ref

cat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/km.ref.tac.dat

• km.hb.tac.dat: High binding TAC.  Average TAC computed across all voxels with labels specified by --km-ref

cat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/km.hb.tac.dat

• gtm.nii.gz: PVCed TACs for each ROI in seg.  This is the input for ROI-based MRTM1 and MRTM2 (not viewable in freeview)
• gtm.stats.dat: Text file with some summary statistics.  Not used for further processing
• mgx.subctxgm.nii.gz: PVC'd TACs (using extended Muller-Gartner) for subcortical areas
• mgx.ctxgm.nii.gz:  PVC'd TACs (using extended Muller-Gartner) for cortical areas

Vis:

cd $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/
freeview \
  --volume ./pet/pet.mc.nii.gz \
  --volume ./pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mgx.subctxgm.nii.gz \
  --volume ./pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mgx.ctxgm.nii.gz

• Use shift+middle mouse click to change contrast
• View --> show time course to vis the TAC at a particular voxel

Creating time.dat file

This is an ASCII text file with the acquisition time (in seconds) of each time point in the time activity curve (TAC).  If your data is in BIDS format, this information is available in the json file associates with the PET data (sub-PS19_ses-baselinebrain_rec-DynTOF_pet.json in the field FrameTimesStart

ls -l $INPUT_DATADIR/ds004230/sub-PS19/ses-baselinebrain/pet
cat $INPUT_DATADIR/ds004230/sub-PS19/ses-baselinebrain/pet/sub-PS19_ses-baselinebrain_rec-DynTOF_pet.json

Extract values to a text file (you can also do this using a text editor like gedit)

cat $INPUT_DATADIR/ds004230/sub-PS19/ses-baselinebrain/pet/sub-PS19_ses-baselinebrain_rec-DynTOF_pet.json \
  | grep -A 33 FrameTimesStart \
  | grep -v FrameTimesStart \
  | sed 's/^ *//' \
  | sed 's/,//' > $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/time.dat

ROI-based MRTM1

mri_glmfit \
  --y $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/km.hb.tac.nii.gz \
  --mrtm1 $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/km.ref.tac.dat \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/time.dat \
  --o $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1 \
  --nii.gz

Where:

• --y input data: TACs for each ROI (output of mri_gtmpvc, we could have also used gtm.nii.gz to compute for each ROI)
• --mrtm1 ...: perform mrtm1 using km.ref.tac.dat as the reference TAC
• $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/time.dat: timing file (created above)
• --o: output folder
• --nii.gz: save output in nii.gz format

View outputs:

ls -l $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1

Outputs:

• bp.nii.gz: Binding potential (for each ROI, here we only have 1 since the input was km.hb.tac.nii.gz)
• k2/gamma.nii.gz: the k2 parameter
• k2a/gamma.nii.gz: the k2a parameter
• R1/gamma.nii.gz: The R1 parameter
• k2prime.dat: The k2' parameter, used as input to MRTM2

View dimensions of outputs:

find \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1 \
  -maxdepth 1 \
  -name '*.nii.gz' \
  -exec sh -c 'echo {}; mri_info --dim {}; echo' \;

None of these files are viewable with FreeView, but they can be converted to text:

mri_convert \
  --ascii \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1/bp.nii.gz \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1/bp.txt
cat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1/bp.txt

ROI-based MRTM2

Let's set an environment variable for k2prime

export K2PRIME=`cat $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm1/k2prime.dat`

Test:

echo $K2PRIME

Now run MRTM2 on all ROIs

mri_glmfit \
  --y $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/gtm.nii.gz \
  --mrtm2 $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/km.ref.tac.dat \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/time.dat \
  $K2PRIME \
  --o $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm2 \
  --nii.gz

View outputs:

ls -l $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm2

View dimensions of outputs:

find \
  $SUBJECTS_DIR/sub-PS19-ses-baselinebrain/pet/gtmpvc.kmref-8-47.kmhb-11-12-13-50-51-52-output/mrtm2 \
  -maxdepth 1 \
  -name '*.nii.gz' \
  -exec sh -c 'echo {}; mri_info --dim {}; echo' \;

Next Steps

Spaces and Group-level Analysis

FreeSurfer outputs are in subject space

Verify:

freeview \
  --volume $INPUT_DATADIR/ds004230/derivatives/freesurfer/sub-PS19/mri/aseg.mgz \
  --volume $INPUT_DATADIR/ds004230/derivatives/freesurfer/sub-PS20/mri/aseg.mgz

But recon-all generates transforms to a common space   - ./mri/transforms/talairach.lta: Linear transform to talairach space   - ./mri/transforms/talairach.m3z: non-linear transform to talairach space

View in talairach space:

freeview \
  --volume $INPUT_DATADIR/ds004230/derivatives/freesurfer/sub-PS19/mri/aseg.mgz:reg=$INPUT_DATADIR/ds004230/derivatives/freesurfer/sub-PS19/mri/transforms/talairach.lta \
  --volume $INPUT_DATADIR/ds004230/derivatives/freesurfer/sub-PS20/mri/aseg.mgz:reg=$INPUT_DATADIR/ds004230/derivatives/freesurfer/sub-PS20/mri/transforms/talairach.lta

See also:

• See sections 7 and 8 on the PetSurfer wiki page
• mri_vol2vol
• mri_vol2surf
• Teach yourself FreeSurfer
  • An Overview of Registration Methods
  • Surface-based analysis: Intersubject Registration & Smoothing
  • Group Analysis
  • Multiple Comparisons
• mris_preproc

Minimamist Defacing (MiDeFace)

See the MiDeFace wiki page

Blood processing

TBD

More Resources

• PetSurfer on the FreeSurfer Wiki
• PetSurfer tutorial
• PetSurfer Overview by Doug Greve
• Teach yourself FreeSurfer
• Pet Processing pipelines
• recon-all dev table
• FreeSurfer Support via Mailinglist
• Example Dataset