Commit 4806dc1f authored by Matteo Bastiani's avatar Matteo Bastiani
Browse files

second data release processing pipeline

parent 2afbf887
......@@ -186,7 +186,7 @@
same "printed page" as the copyright notice for easier
identification within third-party archives.
Copyright 2018 Matteo Bastiani
Copyright 2017 University of Oxford
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
......
-----------------------------------------------
dHCP neonatal dMRI data processing pipeline
March, 2018
-----------------------------------------------
dHCP neonatal dMRI data processing pipeline
April, 2018
V 0.0.2: Processing pipeline used for the 2nd public release
-----------------------------------------------
Automated pipeline to consistently analyse neonatal dMRI data from the developing Human Connectome Project (dHCP).
If you use the pipeline in your work, please cite the following article:
Bastiani, M., Andersson, J.L.R., Cordero-Grande, L., Murgasova, M.,
Hutter, J., Price, A.N., Makropoulos, A., Fitzgibbon, S.P., Hughes,
E., Rueckert, D., Victor, S., Rutherford, M., Edwards, A.D., Smith,
S., Tournier, J.-D., Hajnal, J.V., Jbabdi, S., Sotiropoulos,
S.N. (2019). Automated processing pipeline for neonatal diffusion MRI
in the developing Human Connectome Project. NeuroImage, 185, 750-763.
Installation
------------
The pipeline consists of several bash scripts that do not require any installation.
Once it has been downloaded and unpacked, the first thing to do is fill the necessary paths into the file:
setup.sh
After that, source the script from the terminal in the following way:
. setup.sh
V 0.0.1: Pipeline reflecting the first data release processing
-----------------------------------------------
Comprehensive and automated pipeline to consistently analyse neonatal dMRI data from the developing Human Connectome Project (dHCP).
If you use the pipeline in your work, please cite the following article:
Bastiani, M., Andersson, J., Cordero-Grande, L., Murgasova, M., Hutter, J., Price, A.N., Makropoulos, A., Fitzgibbon, S.P., Hughes, E., Rueckert, D., Victor, S., Rutherford, M., Edwards, A.D., Smith, S., Tournier, J.-D., Hajnal, J.V., Jbabdi, S., Sotiropoulos, S.N. (2018). Automated processing pipeline for neonatal diffusion MRI in the developing Human Connectome Project. NeuroImage.
The script needs to be run before launching the processing jobs.
Installation
Dependencies
------------
The pipeline consists of several bash scripts.
Once it has been downloaded, the first thing to do is fill the correct paths into:
dHCP_neo_dMRI_setup.sh
The script needs to be run before launching the processing jobs.
The dMRI neonatal pipeline mainly relies on:
- FSL 6.0.1 (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki)
Additional dependencies are:
- The dHCP structural pipeline (https://github.com/BioMedIA/dhcp-structural-pipeline)
- ANTs (http://stnava.github.io/ANTs/): non-linear registration to template space
- Convert3D (http://www.itksnap.org/pmwiki/pmwiki.php?n=Convert3D.Documentation): convert ANTs to FNIRT warp fields
- IRTK (https://www.doc.ic.ac.uk/~dr/software/usage.html): super resolution
-------------------------------------------
Examples
-------------------------------------------
To launch the pipeline for a single subject, use the following command:
${scriptsFolder}/dHCP_neo_dMRI_setJobs.sh ${rawDataFolder}/sub-${cid} ses-${no} ${rawDataFile} ${cid} ${scriptsFolder}/dhcp300_f.txt ${outFolder} ${age} ${birth} ${n_sessions}
-------------------------------------------
dHCP data
-------------------------------------------
To launch the pipeline, use:
${scriptsFolder}/dHCP_neo_dMRI_launchPipeline.sh participants.tsv ${outFolder}
This command will submit all the necessary scripts to process the raw dMRI data. The necessary inputs are:
${rawDataFolder}/sub-${cid}: Path to raw data
ses-${no}: Session number
${rawDataFile}: Raw data file name
${cid}: Connectome ID
${scriptsFolder}/dhcp300_f.txt: Protocol file
This command will submit all the necessary jobs to process the dHCP neonatal dMRI datasets.
The two inputs are:
participants.tsv: Subject list
${outFolder}: Output folder
${age}: Age at scan (in rounded weeks)
${birth}: Age at birth (in rounded weeks)
${n_sessions}: Total number of scanning sessions for the same subject
-------------------------------------------
Non-dHCP data
-------------------------------------------
To convert a non-dHCP dataset such that it can be used with the pipeline, use the command:
${scriptsFolder}/utils/getProtocol
By typing the command in a terminal, the necessary inputs will be shown.
The script will generate a raw data file and a protocol file that can be used with the pipeline.
-------------------------------------------
Non-dHCP data
-------------------------------------------
To convert a locally acquired dataset such that it can be used with the pipeline, use the command:
${scriptsFolder}/utils/getProtocol
By typing the command in a terminal, the necessary inputs will be shown.
The script will generate a raw data file and a protocol file that can be used with the pipeline.
Directory structure
-------------------
The pipeline expects the following structure for the input dMRI data:
/path/to/data
/subject1
/session-1
/DWI
/data.nii.gz
/session-2
/DWI
/data.nii.gz
/subject2
/session-1
/DWI
/data.nii.gz
.
.
.
/subjectN
/session-1
/DWI
/data.nii.gz
The getProtocol command can be used to obtain each individual's data.nii.gz 4D nifti volume. All of them will need to be placed in the correct subject/session/DWI folder. This directory structure accounts for the fact that data from a single subject can be acquired in multiple sessions.
Examples
--------
To launch the pipeline for a single subject, use the following command:
${scriptsFolder}/dHCP_neo_dMRI_setJobs.sh ${rawDataFolder}/sub-${cid} \
session-${no} \
${rawDataFile} \
${cid} \
${scriptsFolder}/protocol.txt
${slspec} \
${outFolder} \
${age} \
${birth} \
${subjT2} \
${subjSeg} \
${srFlag} \
${gpuFlag}
This command will submit all the necessary scripts to process the raw dMRI data. Typing the command in the terminal without any input will show the user guide.
The necessary inputs are:
${rawDataFolder}/sub-${cid}: Path to raw data
session-${no}: Session folder
${rawDataFile}: Raw data file name
${cid}: Connectome ID
${scriptsFolder}/protocol.txt: Protocol file
${slspec}: eddy slspec file (0 if not available)
${outFolder}: Output folder
${age}: Age at scan (in weeks, rounded)
${birth}: Age at birth (in weeks, rounded)
${subjT2}: Subject's anatomical T2-weighted volume
${subjSeg}: Subject's tissue segmentation
${srFlag}: super resolution flag (0=do not use super resolution, 1=use
super resolution)
${gpuFlag}: if you have an NVIDIA GPU, this will significantly speed processing time
and allow to use all the eddy features (i.e., slice-to-volume
correction, motion-by-susceptibility-induced distortions)
#!/bin/bash
set -e
echo -e "\n START: getMasks"
unset POSIXLY_CORRECT
if [ "$2" == "" ];then
echo ""
echo "usage: $0 <SubjFolder> <TissueLabels>"
echo " Tissue and brain masks extraction script"
echo " SubjFolder: Path to the subject processing folder"
echo " TissueLabels: Volume containing results of tissue segmentations (assumes 1=CSF, 2=GM, 3=WM, 5=CSF)"
echo ""
exit 1
fi
subjOutFolder=$1 # Path to the subject folder
segVolume=$2 # Segmented volume (assumes 1=CSF, 2=GM, 3=WM, 5=CSF)
anatFolder=${subjOutFolder}/T2w
mkdir -p ${anatFolder}/segmentation
#============================================================================
# Get single tissue masks in T2w space
#============================================================================
${FSLDIR}/bin/imcp ${segVolume} ${anatFolder}/segmentation/tissue_labels.nii # Copy segmented volume to processing folder
${FSLDIR}/bin/fslmaths ${anatFolder}/segmentation/tissue_labels.nii -thr 1 -uthr 1 -bin ${anatFolder}/segmentation/csf_mask # CSF
${FSLDIR}/bin/fslmaths ${anatFolder}/segmentation/tissue_labels.nii -thr 5 -uthr 5 -bin -add ${anatFolder}/segmentation/csf_mask -bin ${anatFolder}/segmentation/csf_mask # Adding ventricles
${FSLDIR}/bin/fslmaths ${anatFolder}/segmentation/tissue_labels.nii -thr 2 -uthr 2 -bin ${anatFolder}/segmentation/gm_mask # GM
${FSLDIR}/bin/fslmaths ${anatFolder}/segmentation/tissue_labels.nii -thr 3 -uthr 3 -bin ${anatFolder}/segmentation/wm_mask # WM
${FSLDIR}/bin/fslmaths ${anatFolder}/segmentation/wm_mask -ero -sub ${anatFolder}/segmentation/wm_mask -abs ${anatFolder}/segmentation/wm_mask_edges.nii.gz # WM edges
#============================================================================
# Create brain mask
#============================================================================
${FSLDIR}/bin/fslmaths ${anatFolder}/segmentation/tissue_labels.nii -thr 0 -bin ${anatFolder}/segmentation/brain_mask
echo -e "\n END: getMasks"
......@@ -4,14 +4,13 @@
set -e
echo -e "\n START: Importing files"
# scriptsFolder=/home/fs0/matteob/scripts/dHCP
if [ "$6" == "" ];then
echo ""
echo "usage: $0 <Data folder> <Data filename> <Output preprocessing folder> <Acquisition protocol> <Number of volumes> <Number of b0s>"
echo "usage: $0 <Data folder> <Data filename> <Output folder> <Acquisition protocol> <Number of volumes> <Number of b0s>"
echo " Data folder: Path to the data file"
echo " Data filename: File name for the raw data"
echo " Output preprocessing folder: Path to the output pre-processing folder"
echo " Output folder: Path to the subject output folder"
echo " Acqusition protocol: Text file with acquisition parameters"
echo " Number of volumes: Number of acquired volumes"
echo " Number of b0s: Number of b0s from each phase encoding block to use when estimating the fieldmap"
......@@ -22,11 +21,15 @@ fi
dataFolder=$1 # Path to data file
dataFile=$2 # Data file name
prepFolder=$3 # Output folder
outFolder=$3 # Subject output folder
acqProt=$4 # Acquisition protcol
nVols=$5 # Number of acquired volumes
noB0s=$6 # Number of B0 volumes for each PE direction used to estimate distortions with TOPUP
prepFolder=${outFolder}/PreProcessed
rawFolder=${outFolder}/raw
mkdir -p ${rawFolder}/tmpData
#============================================================================
# Separate b0 volumes according to their PE direction and obtain
......@@ -47,9 +50,11 @@ no_PA=0
idxb400=0
idxb1000=0
idxb2600=0
bshells=(100000)
n_b=0
i=0
echo -n > ${prepFolder}/eddy/eddyIndex.txt
echo -n > ${rawFolder}/tmpData/eddyIndex.txt
while read line; do
bvec_x[i]=`echo $line | awk {'print $1'}`
bvec_y[i]=`echo $line | awk {'print $2'}`
......@@ -57,55 +62,67 @@ while read line; do
bval[i]=`echo $line | awk {'print $4'}`
pedir[i]=`echo $line | awk {'print $5'}`
rotime[i]=`echo $line | awk {'print $6'}`
echo -n "${pedir[i]} " >> ${prepFolder}/eddy/eddyIndex.txt
if [ ${bval[i]} -eq 0 ]; then
if [ ${pedir[i]} -eq 1 ]; then #LR
LRvols[idxLR]=$i
idxLR=$(($idxLR + 1))
elif [ ${pedir[i]} -eq 2 ]; then #RL
RLvols[idxRL]=$i
idxRL=$(($idxRL + 1))
elif [ ${pedir[i]} -eq 3 ]; then #AP
APvols[idxAP]=$i
idxAP=$(($idxAP + 1))
elif [ ${pedir[i]} -eq 4 ]; then #PA
PAvols[idxPA]=$i
idxPA=$(($idxPA + 1))
fi
else
if [ ${pedir[i]} -eq 1 ]; then #LR
no_LR=$(($no_LR + 1))
elif [ ${pedir[i]} -eq 2 ]; then #RL
no_RL=$(($no_RL + 1))
elif [ ${pedir[i]} -eq 3 ]; then #AP
no_AP=$(($no_AP + 1))
elif [ ${pedir[i]} -eq 4 ]; then #PA
no_PA=$(($no_PA + 1))
fi
if [ ${bval[i]} -eq 400 ]; then
idxb400=$(($idxb400+1))
elif [ ${bval[i]} -eq 1000 ]; then
idxb1000=$(($idxb1000+1))
elif [ ${bval[i]} -eq 2600 ]; then
idxb2600=$(($idxb2600+1))
fi
echo -n "${pedir[i]} " >> ${rawFolder}/tmpData/eddyIndex.txt
if [ ${bval[i]} -lt 100 ]; then #b0
first_b0=${i}
if [ ${pedir[i]} -eq 1 ]; then #LR
LRvols[idxLR]=$i
idxLR=$(($idxLR + 1))
elif [ ${pedir[i]} -eq 2 ]; then #RL
RLvols[idxRL]=$i
idxRL=$(($idxRL + 1))
elif [ ${pedir[i]} -eq 3 ]; then #AP
APvols[idxAP]=$i
idxAP=$(($idxAP + 1))
elif [ ${pedir[i]} -eq 4 ]; then #PA
PAvols[idxPA]=$i
idxPA=$(($idxPA + 1))
fi
else #dw
if [ ${pedir[i]} -eq 1 ]; then #LR
no_LR=$(($no_LR + 1))
elif [ ${pedir[i]} -eq 2 ]; then #RL
no_RL=$(($no_RL + 1))
elif [ ${pedir[i]} -eq 3 ]; then #AP
no_AP=$(($no_AP + 1))
elif [ ${pedir[i]} -eq 4 ]; then #PA
no_PA=$(($no_PA + 1))
fi
fi
# Identify unique shells
b_flag=0
for ub in "${bshells[@]}"; do
j=`echo ${bval[i]} | awk -F"E" 'BEGIN{OFMT="%10.10f"} {print $1 * (10 ^ $2)}' ` # Round bval
j=${j/.*}
b_diff=`echo "${j} - ${ub}" | bc | awk ' { if($1>=0) { print $1} else {print $1*-1 }}'`
if [ ${b_diff} -lt 100 ]; then
b_flag=1
break
fi
done
if [ "${b_flag}" == 0 ]; then
bshells[n_b]=${bval[i]}
n_b=$((${n_b} + 1))
fi
i=$(($i + 1))
if [ $i -eq ${nVols} ]; then
break
break
fi
done < ${acqProt}
bshells=($(echo "${bshells[@]}" | tr ' ' '\n' | sort -n -u | tr '\n' ' '))
echo "${bshells[@]}" > ${rawFolder}/shells
echo "Found ${n_b} shells: ${bshells[@]} s/mm^2"
un_pedirs=(`echo "${pedir[@]}" | tr ' ' '\n' | sort -n -u | tr '\n' ' '`)
echo "${un_pedirs[@]}" > ${rawFolder}/pedirs
#============================================================================
# Write json file for QC
#============================================================================
echo "{" > ${prepFolder}/dataImport.json
echo " \"no_B400_vols\": $idxb400," >> ${prepFolder}/dataImport.json
echo " \"no_B1000_vols\": $idxb1000," >> ${prepFolder}/dataImport.json
echo " \"no_B2600_vols\": $idxb2600," >> ${prepFolder}/dataImport.json
echo " \"no_LR_vols\": $no_LR," >> ${prepFolder}/dataImport.json
echo " \"no_RL_vols\": $no_RL," >> ${prepFolder}/dataImport.json
echo " \"no_AP_vols\": $no_AP," >> ${prepFolder}/dataImport.json
......@@ -117,10 +134,10 @@ echo "}" >> ${prepFolder}/dataImport.json
#============================================================================
# Write bvecs, bvals and eddy acquisition parameters file.
#============================================================================
echo "${bvec_x[@]}" > ${prepFolder}/tmpData/bvecs
echo "${bvec_y[@]}" >> ${prepFolder}/tmpData/bvecs
echo "${bvec_z[@]}" >> ${prepFolder}/tmpData/bvecs
echo "${bval[@]}" > ${prepFolder}/bvals
echo "${bvec_x[@]}" > ${rawFolder}/tmpData/orig_bvecs
echo "${bvec_y[@]}" >> ${rawFolder}/tmpData/orig_bvecs
echo "${bvec_z[@]}" >> ${rawFolder}/tmpData/orig_bvecs
echo "${bval[@]}" > ${rawFolder}/tmpData/bvals
if [ $LRvols -ne -1 ]; then
echo -1 0 0 ${rotime[${LRvols[0]}]} > ${prepFolder}/eddy/acqparamsUnwarp.txt
......@@ -143,95 +160,41 @@ else
echo 0 1 0 0.05 >> ${prepFolder}/eddy/acqparamsUnwarp.txt
fi
#============================================================================
# Reorient raw data and bvecs to standard orientations
# Reorient raw data and bvecs to to match the approximate orientation
# of the standard template images (MNI152)
#============================================================================
${FSLDIR}/bin/fslreorient2std ${dataFolder}/${dataFile} ${prepFolder}/tmpData/data
${FSLDIR}/bin/fslreorient2std ${dataFolder}/${dataFile} > ${prepFolder}/tmpData/raw2std.mat
${FSLDIR}/bin/fslreorient2std ${dataFolder}/${dataFile} ${rawFolder}/tmpData/data
${FSLDIR}/bin/fslreorient2std ${dataFolder}/${dataFile} > ${rawFolder}/tmpData/raw2std.mat
echo 1 0 0 0 > ${prepFolder}/tmpData/flipZ.mat
echo 0 1 0 0 >> ${prepFolder}/tmpData/flipZ.mat
echo 0 0 -1 0 >> ${prepFolder}/tmpData/flipZ.mat
echo 0 0 0 1 >> ${prepFolder}/tmpData/flipZ.mat
${scriptsFolder}/utils/rotateBvecs.sh ${rawFolder}/tmpData/orig_bvecs ${rawFolder}/tmpData/raw2std.mat ${rawFolder}/tmpData/bvecs
${scriptsFolder}/utils/rotateBvecs.sh ${prepFolder}/tmpData/bvecs ${prepFolder}/tmpData/raw2std.mat ${prepFolder}/tmpData/bvecs2
${scriptsFolder}/utils/rotateBvecs.sh ${prepFolder}/tmpData/bvecs2 ${prepFolder}/tmpData/flipZ.mat ${prepFolder}/bvecs
#============================================================================
# Check that in-plane matrix size is a multiple of 2 (for TOPUP)
# If more than 1 PE direction has been acquired, Identify best b0s for each
# PE direction; otherwise, set the first b0 as the reference volume.
#============================================================================
dimt1=`${FSLDIR}/bin/fslval ${prepFolder}/tmpData/data dim1`
c1=$(($dimt1%2))
if [ $c1 -ne 0 ]; then
${FSLDIR}/bin/fslroi ${prepFolder}/tmpData/data ${prepFolder}/tmpData/tmp 0 1 0 -1 0 -1 0 -1
${FSLDIR}/bin/fslmerge -x ${prepFolder}/tmpData/data ${prepFolder}/tmpData/data ${prepFolder}/tmpData/tmp
fi
dimt2=`${FSLDIR}/bin/fslval ${prepFolder}/tmpData/data dim2`
c2=$(($dimt2%2))
if [ $c2 -ne 0 ]; then
${FSLDIR}/bin/fslroi ${prepFolder}/tmpData/data ${prepFolder}/tmpData/tmp 0 -1 0 1 0 -1 0 -1
${FSLDIR}/bin/fslmerge -y ${prepFolder}/tmpData/data ${prepFolder}/tmpData/data ${prepFolder}/tmpData/tmp
unique_pedirs=(`cat ${rawFolder}/pedirs`)
if [ `echo ${#unique_pedirs[@]}` -gt 1 ]; then
echo "More than 1 phase encoding direction detected. Selecting best b0 volumes for topup"
${scriptsFolder}/utils/pickBestB0s ${rawFolder}/tmpData/data ${rawFolder}/tmpData/bvals ${rawFolder}/tmpData/eddyIndex.txt ${rotime[0]} ${noB0s} ${prepFolder}/topup
else
echo "1 phase encoding direction detected. Setting first b0 as reference volume"
echo "${first_b0}" > ${prepFolder}/topup/ref_b0.txt
fi
#============================================================================
# Identify best b0s for every acquired PE direction
# Sort raw data based on the selected reference volume
#============================================================================
echo "${LRvols[@]}" > ${prepFolder}/topup/b0Indices
echo "${RLvols[@]}" >> ${prepFolder}/topup/b0Indices
echo "${APvols[@]}" >> ${prepFolder}/topup/b0Indices
echo "${PAvols[@]}" >> ${prepFolder}/topup/b0Indices
if [ $LRvols -ne -1 ]; then
tmp=$(printf ",%s" "${LRvols[@]}")
tmp=${tmp:1}
${FSLDIR}/bin/fslselectvols -i ${prepFolder}/tmpData/data -o ${prepFolder}/topup/LR_B0s --vols=${tmp}
$scriptsFolder/utils/pickBestB0s.sh ${prepFolder}/topup 0 ${noB0s}
fi
if [ $RLvols -ne -1 ]; then
tmp=$(printf ",%s" "${RLvols[@]}")
tmp=${tmp:1}
${FSLDIR}/bin/fslselectvols -i ${prepFolder}/tmpData/data -o ${prepFolder}/topup/RL_B0s --vols=${tmp}
$scriptsFolder/utils/pickBestB0s.sh ${prepFolder}/topup 1 ${noB0s}
fi
if [ $APvols -ne -1 ]; then
tmp=$(printf ",%s" "${APvols[@]}")
tmp=${tmp:1}
${FSLDIR}/bin/fslselectvols -i ${prepFolder}/tmpData/data -o ${prepFolder}/topup/AP_B0s --vols=${tmp}
$scriptsFolder/utils/pickBestB0s.sh ${prepFolder}/topup 2 ${noB0s}
fi
if [ $PAvols -ne -1 ]; then
tmp=$(printf ",%s" "${PAvols[@]}")
tmp=${tmp:1}
${FSLDIR}/bin/fslselectvols -i ${prepFolder}/tmpData/data -o ${prepFolder}/topup/PA_B0s --vols=${tmp}
$scriptsFolder/utils/pickBestB0s.sh ${prepFolder}/topup 3 ${noB0s}
fi
${scriptsFolder}/utils/sortData ${rawFolder}/tmpData/data `cat ${prepFolder}/topup/ref_b0.txt` ${rawFolder} ${rawFolder}/tmpData/bvals ${rawFolder}/tmpData/bvecs ${rawFolder}/tmpData/eddyIndex.txt
#============================================================================
# Sort b0s based on their scores, merge them in a file and write acqparams
# for topup and reference scan for eddy.
# Clean unnecessary files
#============================================================================
bestPE=(`cat ${prepFolder}/topup/idxBestB0s.txt | sort -k 3 -n | head -n 4 | awk '{print $1}'`)
bestVol=(`cat ${prepFolder}/topup/idxBestB0s.txt | sort -k 3 -n | head -n 4 | awk '{print $2}'`)
acqpPE=("-1 0 0" "1 0 0" "0 -1 0" "0 1 0")
PE=("LR" "RL" "AP" "PA")
count=0
echo -n > ${prepFolder}/topup/acqparams.txt
for i in "${bestPE[@]}" ; do
dimt4=`${FSLDIR}/bin/fslval ${prepFolder}/topup/"${PE[$i]}"_B0s.nii.gz dim4`
for j in $(seq 0 $((${dimt4}-1))) ; do
echo "${acqpPE[$i]}" "${rotime[${bestVol[${count}]}]}" >> ${prepFolder}/topup/acqparams.txt
done
PE_B0s[$count]=${prepFolder}/topup/"${PE[$i]}"_B0s.nii.gz
count=$(($count+1))
done
echo ${PE_B0s[@]}
${FSLDIR}/bin/fslmerge -t ${prepFolder}/topup/phase ${PE_B0s[@]}
echo "${bestVol[0]}" > ${prepFolder}/eddy/ref_scan.txt
rm -rf ${rawFolder}/tmpData
echo -e "\n END: Importing files"
......
#!/bin/bash
echo "\n START: dHCP neonatal dMRI data processing pipeline"
if [ "${2}" == "" ];then
echo "The script will read dHCP subject info and, if data is there, launch the processing steps"
echo ""
echo "usage: dHCP_neo_dMRI_launchPipeline.sh <subject list> <output folder>"
echo ""
echo " subject list: text file containing participant_id, sex and age at birth (w GA)"
echo " output folder: folder where results will be stored"
echo ""
echo ""
fi
subjList=$1
outFolder=$2
mkdir -p ${outFolder}
# Read the connectome IDs
sids=(`cat ${subjList} | sed "1 d" | cut -f 1 | grep -v "^$"`)
# Main loop through subjects
for s in ${sids[@]}; do
sex=`cat ${subjList} | grep ${s} | cut -f 2`
birth=`cat ${subjList} | grep ${s} | cut -f 3` # Age at birth (weeks)
sessions=(`cat ${reconFolder}/sub-${s}/sub-${s}_sessions.tsv | sed "1 d" | cut -f 1 | grep -v "^$"`) # Some subjects are acquired over multiple sessions
n_sessions=`echo ${#sessions[@]}`
if [ ${n_sessions} -gt 0 ]; then
for ses in ${sessions[@]}; do
date=`cat ${reconFolder}/sub-${s}/sub-${s}_sessions.tsv | grep ${ses} | cut -f 2`
age=`cat ${reconFolder}/sub-${s}/sub-${s}_sessions.tsv | grep ${ses} | cut -f 3` # Age at scan (weeks)
# Round birth and scan ages
age=`awk -v v="${age}" 'BEGIN{printf "%.0f", v}'`
birth=`awk -v v="${birth}" 'BEGIN{printf "%.0f", v}'`
# Subject specific variables
data=${reconFolder}/sub-${s}/ses-${ses}/DWI/sub-${s}_ses-${ses}_DWI_MB0_AnZfZfAdGhAb.nii
t2=${structFolder}/sub-${s}/ses-${ses}/anat/sub-${s}_ses-${ses}_T2w_restore.nii.gz
seg=${structFolder}/sub-${s}/ses-${ses}/anat/sub-${s}_ses-${ses}_drawem_tissue_labels.nii.gz
if [ -e ${data} ]; then # Check that data has been acquired
#============================================================================
# Check for scan completeness
#============================================================================
dimt4=`${FSLDIR}/bin/fslval ${data} dim4`
complete_check=${dimt4}
usable_check=1
if [ ${dimt4} -lt 34 ]; then
echo "WARNING: The dataset is unusable as it does not contain enough b0 volumes"
echo "${s} ses-${ses}" >> ${outFolder}/unusable.txt
usable_check=0
elif [ ${dimt4} -lt 123 ]; then
echo "WARNING: The dataset is incomplete and does not contain enough b0 pairs for each PE direction"
echo "${s} ses-${ses}" >> ${outFolder}/incomplete.txt
noB0s=1
usable_check=1
fi
#============================================================================
# Store QC information
#============================================================================
subjOutFolder=${outFolder}/${s}/ses-${ses}
if [ -e ${subjOutFolder}/initQC.json ]; then
break
fi
mkdir -p ${subjOutFolder}
echo "{" > ${subjOutFolder}/initQC.json
echo " \"Complete\": ${complete_check}," >> ${subjOutFolder}/initQC.json
echo " \"Usable\": ${usable_check}," >> ${subjOutFolder}/initQC.json
echo " \"nSessions\": ${n_sessions}," >> ${subjOutFolder}/initQC.json
echo " \"birthAge\": ${birth}," >> ${subjOutFolder}/initQC.json
echo " \"scanAge\": ${age}" >> ${subjOutFolder}/initQC.json
echo "}" >> ${subjOutFolder}/initQC.json
if [ -e ${t2} ]; then # Check that structural data has been acquired
#============================================================================
# Set processing jobs
#============================================================================
${scriptsFolder}/dHCP_neo_dMRI_setJobs.sh ${reconFolder}/sub-${s} ses-${ses} sub-${s}_ses-${ses}_DWI_MB0_AnZfZfAdGhAb.nii ${s} \
${scriptsFolder}/dHCP_protocol.txt ${scriptsFolder}/slorder.txt ${outFolder} \
${age} ${birth} ${t2} ${seg} 1 1
echo "${s} ${ses} ${birth} ${age}" >> ${outFolder}/complete.txt
else
echo "WARNING! Missing structural data for subject ${s}"
echo "${s} ses-${ses}" >> ${outFolder}/missingAnat.txt
fi
else
echo "WARNING! Missing dMRI data for subject ${s}"
echo "${s} ses-${ses}" >> ${outFolder}/missingDmri.txt
fi
done
else
echo "WARNING! Missing session IDs for subject ${s}"
echo "${s}" >> ${outFolder}/missingSessions.txt
fi
done
echo "\n END: dHCP neonatal dMRI data processing pipeline"