Computing local gyrification from CIVET outputs

make-gyrification.sh is a post-processing script that will extract vertex-wise local gyrification index from surfaces generated by CIVET.

Below is the script that will generate a gyrification_joblist, which you can then submit to the cluster through qbatch to execute the make-gyrification commands for each subject in parallel. Depending on the version of CIVET you used to process your T1-weighted images and your own directory names, some lines of the code will need to be changed (see comments below).

Modules required (Niagara)

> module load cobralab/2019b # if you did not use CIVET/2.1.1, you will need to unload it from the list and load the version that you used + dependencies 

Working directory structure

make-gyrification.sh should be run from the directory containing your CIVET output folder. As well, the script needs access to the icbm reference left/right hemisphere surfaces used in CIVET. On Niagara, these two surfaces can be found in:

• $QUARANTINE_PATH/CIVET/2.1.0/src/Linux-x86_64/SRC/civet-2.1.0/models/icbm/ (for CIVET 2.1)
• $QUARANTINE_PATH/resources/CIVET/ (for CIVET 2.0)

The reference surfaces will be .obj files with "icbm_avg_mid_sym_mc_{left,right}" in the name - copy these to your local directory containing make-gyrification.sh and your CIVET output folder. The script is not optimized for hires outputs from CIVET (163,842 vertices/hemi), so ensure that you do not copy the _hires surfaces.

You will also a copy of the surfaceratio scripts, found here (on CIC): /data/chamal/projects/alyssaD/PS_Dimensions/derivatives/scripts/surfaceratio

To summarize, the working directory should contain:

• make-gyrification.sh (see below)
• your CIVET output folder
• the icbm average surfaces for the version of CIVET used
• a surfaceratio/ folder

Running make-gyrification.sh

Make the script executable using the command chmod +x ./make-gyrification.sh, then run it to generate your joblist.

# make-gyrification.sh

# In this example, 2.1.0_ref_surfaces/ is a local folder containing the two reference surfaces. Adjust this for your own use as needed.

for left in output/*/surfaces/*_mid_surface_left_81920.obj; do
	echo "./surfaceratio/surfaceratio ${left} $(dirname $left)/gyrification_left.txt 20 && \
	surface-resample 2.1.0_ref_surfaces/icbm_avg_mid_sym_mc_left.obj \
	${left} \
	$(dirname ${left})/gyrification_left.txt \
	$(dirname $(dirname ${left}))/transforms/surfreg/*_left_surfmap.sm \
	$(dirname ${left})/gyrification_rsl_left.txt" >> gyrification_joblist
done

for right in output/*/surfaces/*_mid_surface_right_81920.obj; do
	echo "./surfaceratio/surfaceratio ${right} $(dirname $right)/gyrification_right.txt 20 && \
	surface-resample 2.1.0_ref_surfaces/icbm_avg_mid_sym_mc_right.obj \
	${right} \
	$(dirname ${right})/gyrification_right.txt \
	$(dirname $(dirname ${right}))/transforms/surfreg/*_right_surfmap.sm \
	$(dirname ${right})/gyrification_rsl_right.txt" >> gyrification_joblist
done 

Once you have your gyrification_joblist, you can submit it to Niagara:

>  qbatch -w 00:15:00 gyrification_joblist # For reference, for 153 subjects this took only a few minutes. Unless your dataset is really huge, the default 15 min should be sufficient.

Outputs

Once your jobs are finished, you should have 4 new files in each of output/*/surfaces/: gyrification_left.txt, gyrification_right.txt, gyrification_rsl_left.txt, gyrification_rsl_right.txt. These contain vertex-wise gyrification measurements, where rsl denotes measurements resampled on the surface model (what we're usually interested in).

Computing regional averages

If you're interested in ROIs rather than vertices, run the following Bash script in the directory containing your CIVET output folder to compute an average local gyrification for each ROI in the DKT or AAL atlases from vertex-wise values from make-gyrification.sh.

Requirements:

• get_regional_gyrif.py
• local copies of the left/right hemisphere atlas files (DKT/AAL) containing a region number for each vertex, also found in $QUARANTINE_PATH/resources/CIVET/ (use the _short versions for 40,962 vertices/hemi)
• local copy of the atlas region number-to-name mapping, found in $QUARANTINE_PATH/CIVET/2.1.0/install/CIVET-2.1.0/models/icbm/DKT/DKTatlas40.labels for the DKT atlas, and $QUARANTINE_PATH/resources/CIVET/CIVET_2.0_AAL_labels.txt for the AAL atlas
• cobralab/2019b module

#!/bin/bash

# make_regional_gyrif.sh
# Desc: Runs get_regional_gyrif.py on ROI-labelled gyrification files of each subject, per hemisphere. Output files are moved into subjects' surfaces/.

## Add column to vertex-wise gyrification text files (from make-gyrification.sh) containing DKT region label associated with each vertex
# Modify below variables for your data path/selected atlas
atlas_left=/path/to/my_civet_directory/DKT/CIVET_2.1.0_dkt_left_short.txt
atlas_right=/path/to/my_civet_directory/DKT/CIVET_2.1.0_dkt_right_short.txt
atlas_names=/path/to/my_civet_directory/DKT/DKTatlas40.labels
subject_dir=/path/to/my_civet_directory/output
feature_left=DKT_gyrification_rsl_left # feature name
feature_right=DKT_gyrification_rsl_right

# The DKT and AAL region number-to-name files are formatted differently. The below code cleans up the AAL region names file to match the DKT file if AAL is used
if echo "$atlas_names" | grep -q "AAL"; then
  is_aal=true
  cat $atlas_names | tr -s "[:blank:]" " " | cut -d" " -f 1,2 | tr -d "'" | tail -n +2 > temp_atlas_names.txt
  atlas_names=temp_atlas_names.txt
fi

# Split atlases names file into left/right hemis
head -n 32 ${atlas_names} > DKT_names_L.txt
tail -n 32 ${atlas_names} > DKT_names_R.txt
atlas_names_L=DKT_names_L.txt
atlas_names_R=DKT_names_R.txt

for subject in $subject_dir/*/; do
  subID=$(basename $subject)

  paste $atlas_left ${subject}surfaces/gyrification_rsl_left.txt | column -o $'\t' -t > ${subject}surfaces/${subID}_${feature_left}.txt
  paste $atlas_right ${subject}surfaces/gyrification_rsl_right.txt | column -o $'\t' -t > ${subject}surfaces/${subID}_${feature_right}.txt
done


## Run get_regional_gyrif.py
for subject in $subject_dir/*/; do
  subID=$(basename $subject)

  # left hemi
  lh_gyrif=${subject}surfaces/${subID}_${feature_left}.txt
  lh_filename=${subID}_${feature_left}
  # echo $lh_filename
  python3 get_regional_gyrif.py $lh_gyrif $lh_filename $atlas_names_L

  # right hemi
  rh_gyrif=${subject}surfaces/${subID}_${feature_right}.txt
  rh_filename=${subID}_${feature_right}
  # echo $rh_filename
  python3 get_regional_gyrif.py $rh_gyrif $rh_filename $atlas_names_R

  # modify for your file path
  mv /path/to/my_civet_directory/${lh_filename}.csv ${subject}surfaces/
  mv /path/to/my_civet_directory/${rh_filename}.csv ${subject}surfaces/
done

# delete temporary file generated if AAL atlas is chosen
if [ "$is_aal" = true ]; then
  rm temp_atlas_names.txt
fi

# get_regional_gyrification.py

# Desc: Computes average local gyrification for each region of a specified atlas.
# Output: csv file with first column - region number, second column - region name, third column - ROI average gyrification value.

import sys
import pandas as pd

filepath = sys.argv[1]
filename = sys.argv[2]
atlasnamespath = sys.argv[3]

vert_gyrif = pd.read_csv(file, sep='\t', names=['Label', 'ROI value'])
roi_gyrif = vert_gyrif.groupby('Label')['ROI value'].mean().reset_index()
roi_gyrif = roi_gyrif.round(4)

roi_names = pd.read_csv(file, sep=' ', names=['Label', 'ROI name'])
roi_names_gyrif = pd.merge(roi_names, roi_gyrif, on = 'Label', how = 'right')

roi_names_gyrif.to_csv(filename + '.csv', sep=',', index=False, header=False)