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.