Chapter 14 Automatic Fiber Quantification

After you’ve preprocessed your diffusion weighted data using dtiInit, you are ready to run the AFQ pipeline.

14.1 Parameters

AFQ pipeline needs several MATLAB vectors in order to run. First, you will need to set the path to each participant’s dt6.mat file under the variable sub_dirs, then you will need to set the group information (0 for the control group and 1 for study group) under the variable sub_group. I do this arbitrarily because I don’t often have participant ID beforehand. Setting groups therefore is done randomly and doesn’t not effect analyses. I’ve created a script that will generate the variables, sub_dirs and sub_group.

mkdir -p /work/ashley.ware/data/ACAP/afq
workdir=/work/ashley.ware/data/ACAP/afq
rm ${workdir}/sub_dirs.txt
rm ${workdir}/sub_group.txt
rm ${workdir}/control.csv
rm ${workdir}/patient.csv
for i in $(find /work/ashley.ware/data/ACAP/dtiInit/ -type d -name "dt"); do
  subjid=`echo $i | cut -d "/" -f7`
  echo "$i" >> ${workdir}/sub_dirs.txt
  group=$((RANDOM%2))
  if [ $group == 1 ]; then
    echo -n "1, " >> ${workdir}/sub_group.txt
    echo "$subjid, 1" >> ${workdir}/patient.csv
  else
    echo -n "0, " >> ${workdir}/sub_group.txt
    echo "$subjid, 0" >> ${workdir}/control.csv
  fi
done

14.2 Batch Script

First check that your total count of subjects seems correct.

total=`find /work/ashley.ware/data/ACAP/dtiInit/ -type d -name "dt" | wc -l`
echo $total

If everything is correct then you can create your batch file.

Create script.

vi ~/scripts/ACAP/afqIndividual-batch.sh

Copy and paste.

#!/bin/bash
curTime=`date +"%Y%m%d-%H%M%S"`
username=`id -un`
mkdir -p ~/logfiles/ACAP/${curTime}
total=`find /work/${username}/data/ACAP/dtiInit/ -type d -name "dt" | wc -l`
for subjnum in `seq 1 $total`; do
  sbatch \
  -o ~/logfiles/ACAP/${curTime}/output-${subjnum}.txt \
  -e ~/logfiles/ACAP/${curTime}/error-${subjnum}.txt \
  ~/scripts/ACAP/afqIndividual-job.sh \
  ${subjnum}
  sleep 1
done

14.3 Job Script

Create job script.

vi ~/scripts/ACAP/afqIndividual-job.sh

Copy and paste.

#!/bin/bash

#SBATCH --time=03:00:00   # walltime
#SBATCH --ntasks=1  # number of processor cores (i.e. tasks)
#SBATCH --nodes=1   # number of nodes
#SBATCH --mem-per-cpu=16G   # memory per CPU core

# LOAD MODULES, INSERT CODE, AND RUN YOUR PROGRAMS HERE
cd ~/scripts/ACAP/
module load matlab
module load lib/openblas/0.2.20-gnu
module load fsl/6.0.0
matlab -nodisplay -nojvm -nosplash -r "afqIndividual($1)"

14.4 MATLAB Script

Create MATLAB script.

vi ~/scripts/ACAP/afqIndividual.m

Copy and paste.

%%%%%%%
function afqIndividual(x)

% Get home directory:
var = getenv('HOME');

% Add modules to MATLAB. Do not change the order of these programs:
SPM8Path = [var, '/apps/matlab/spm8'];
addpath(genpath(SPM8Path));
vistaPath = [var, '/apps/matlab/vistasoft'];
addpath(genpath(vistaPath));
AFQPath = [var, '/apps/matlab/AFQ'];
addpath(genpath(AFQPath));

% AFQ
sub_dirs = importdata('/work/ashley.ware/data/ACAP/afq/sub_dirs.txt');
load /work/ashley.ware/data/ACAP/afq/sub_group.txt
C = strsplit(char(sub_dirs(x)), '/')
outdir = fullfile(strcat('/work/ashley.ware/data/ACAP/afq/', C(7)));
mkdir(char(outdir));
cd(char(outdir));
if isempty(dir('afq*.mat'))
outfile = strcat('afq-', datestr(now, 'yyyy-mm-dd-HHMM'));
outname = fullfile(outdir, outfile)
afq = AFQ_Create('sub_dirs', sub_dirs, 'sub_group', sub_group, 'showfigs', false);
afq = AFQ_set(afq, 'runsubjects', x);
afq = AFQ_set(afq, 'computenorms', 0);
[afq] = AFQ_run(sub_dirs, sub_group, afq);
save(char(outname), 'afq');

% Corpus Callosum
outfile = strcat('cc-', datestr(now, 'yyyy-mm-dd-HHMM'));
outname = fullfile(outdir, outfile)
afq = AFQ_SegmentCallosum(afq, [false, false]);
save(char(outname), 'afq');

% Export Data
A = afq.vals.fa;
B =  [A{1}(x,:).' A{2}(x,:).' A{3}(x,:).' A{4}(x,:).' A{5}(x,:).' A{6}(x,:).' A{7}(x,:).' A{8}(x,:).' A{9}(x,:).' A{10}(x,:).' A{11}(x,:).' A{12}(x,:).' A{13}(x,:).' A{14}(x,:).' A{15}(x,:).' A{16}(x,:).' A{17}(x,:).' A{18}(x,:).' A{19}(x,:).' A{20}(x,:).' A{21}(x,:).' A{22}(x,:).' A{23}(x,:).' A{24}(x,:).' A{25}(x,:).' A{26}(x,:).' A{27}(x,:).' A{28}(x,:).'];
outfile = strcat('fa.csv');
outname = fullfile(outdir, outfile);
csvwrite(char(outname),B);

A = afq.vals.md;
B =  [A{1}(x,:).' A{2}(x,:).' A{3}(x,:).' A{4}(x,:).' A{5}(x,:).' A{6}(x,:).' A{7}(x,:).' A{8}(x,:).' A{9}(x,:).' A{10}(x,:).' A{11}(x,:).' A{12}(x,:).' A{13}(x,:).' A{14}(x,:).' A{15}(x,:).' A{16}(x,:).' A{17}(x,:).' A{18}(x,:).' A{19}(x,:).' A{20}(x,:).' A{21}(x,:).' A{22}(x,:).' A{23}(x,:).' A{24}(x,:).' A{25}(x,:).' A{26}(x,:).' A{27}(x,:).' A{28}(x,:).'];
outfile = strcat('md.csv');
outname = fullfile(outdir, outfile);
csvwrite(char(outname),B);

A = afq.vals.rd;
B =  [A{1}(x,:).' A{2}(x,:).' A{3}(x,:).' A{4}(x,:).' A{5}(x,:).' A{6}(x,:).' A{7}(x,:).' A{8}(x,:).' A{9}(x,:).' A{10}(x,:).' A{11}(x,:).' A{12}(x,:).' A{13}(x,:).' A{14}(x,:).' A{15}(x,:).' A{16}(x,:).' A{17}(x,:).' A{18}(x,:).' A{19}(x,:).' A{20}(x,:).' A{21}(x,:).' A{22}(x,:).' A{23}(x,:).' A{24}(x,:).' A{25}(x,:).' A{26}(x,:).' A{27}(x,:).' A{28}(x,:).'];
outfile = strcat('rd.csv');
outname = fullfile(outdir, outfile);
csvwrite(char(outname),B);

A = afq.vals.ad;
B =  [A{1}(x,:).' A{2}(x,:).' A{3}(x,:).' A{4}(x,:).' A{5}(x,:).' A{6}(x,:).' A{7}(x,:).' A{8}(x,:).' A{9}(x,:).' A{10}(x,:).' A{11}(x,:).' A{12}(x,:).' A{13}(x,:).' A{14}(x,:).' A{15}(x,:).' A{16}(x,:).' A{17}(x,:).' A{18}(x,:).' A{19}(x,:).' A{20}(x,:).' A{21}(x,:).' A{22}(x,:).' A{23}(x,:).' A{24}(x,:).' A{25}(x,:).' A{26}(x,:).' A{27}(x,:).' A{28}(x,:).'];
outfile = strcat('ad.csv');
outname = fullfile(outdir, outfile);
csvwrite(char(outname),B);
else
sprintf('AFQ is already completed')
end
exit

14.5 Submit Jobs with Batch Script

Submit batch script.

sh ~/scripts/ACAP/afqIndividual-batch.sh

You can check the status of your jobs by checking if it is running or pending.

squeue

You can check the outputs.

tail -n 5 ~/logfiles/ACAP/<date>/output*

14.6 Download Directory / Sync Directory with Arc

Back on your local computer make sure you have the directory already.

mkdir -p /Volumes/bobo/data/ACAP/afq

Then sync the data from the remote computer to your local computer.

rsync -rauv \
ashley.ware@arc.ucalgary.ca:/work/ashley.ware/data/ACAP/afq/ \
/Volumes/bobo/data/ACAP/afq/

14.7 Verify AFQ Analysis

Make directory to save figures:

mkdir -p /Volumes/bobo/data/ACAP/afq/figures/

Open MATLAB and run the following code:

% Get home directory:
var = getenv('HOME');

% Add modules to MATLAB. Do not change the order of these programs:
SPM8Path = [var, '/apps/spm8'];
addpath(genpath(SPM8Path));
vistaPath = [var, '/apps/vistasoft'];
addpath(genpath(vistaPath));
AFQPath = [var, '/apps/afq'];
addpath(genpath(AFQPath));

pathDir = '/Volumes/bobo/data/ACAP/dtiInit/';
subjid = dir(pathDir);

% % Forceps Major and Minor
for i = 1:length(subjid)

  if exist(fullfile([pathDir, subjid(i).name, ...
    '/dt/fibers/MoriGroups_clean_D5_L4.mat']), 'file')
    
    subjid(i).name
    
    % Load Fiber Groups
    fg = dtiReadFibers(fullfile([pathDir, subjid(i).name, ...
    '/dt/fibers/MoriGroups_clean_D5_L4.mat']));
    dt = dtiLoadDt6(fullfile([pathDir, subjid(i).name, ...
      '/dt/dt6.mat']));
    t1 = readFileNifti(fullfile([pathDir, subjid(i).name, ...
      '/dt/bin/wmProb.nii.gz']));
 
    % Create Figure
    AFQ_RenderFibers(fg(9), 'numfibers', 500, ...
      'color', [1 0 0], 'camera', 'sagittal');
    AFQ_RenderFibers(fg(10), 'numfibers', 500, ...
      'color', [1 0 0], 'newfig', '0');
    AFQ_AddImageTo3dPlot(t1, [1, 0, 0]);
 
    % Format Figure
    set(gcf, 'Position', [100, 100, 750, 538]);
    set(gca, 'XTick', [], 'YTick', [], 'ZTick', [], ...
      'xlabel', [], 'ylabel', [], 'zlabel', []);
    set(gca, 'LooseInset', get(gca, 'TightInset'))
 
    % Save Figure
    saveas(gcf, ...
    fullfile([...
    '/Volumes/bobo/data/ACAP/afq/figures/', ...
      subjid(i).name, '.png']));
    close all
  end
end

You’re output images are located at /Volumes/bobo/data/ACAP/afq/figures/ and they should look something like this. Showing the corpus callosum is done to make sure they tensor directions were set correctly. If left and right were flipped, which is common, your corpus callosum tracks would look like a tangled mess. Just an FYI, you can use the MATLAB code above the generate figures of any of the fiber tracks.

sub-ACAP1006_ses-FU6

sub-ACAP1006_ses-FU6