Chapter 6 Diffusion MRI

6.1 DWI pre-processing

#!/bin/bash

## dwi preprocessing
start=`date +%s`

nif_dir=/home/neuro/MRI_GDPH/nifti
proc_dir=/home/neuro/MRI_GDPH/processed

# basename   /home/neuro/MRI_GDPH/nifti/PD-348/anat/PD-348_visit_1_t1.nii.gz _t1.nii.gz
## output: PD-348_visit_1

# dti_list=$(find . -name '*dti.nii*')
# echo "$dti_list" >> /home/neuro/codebank/GDPH/dti_list.txt

#*dti.nii.gz
# file=*dti.nii.gz
## ./PD-343/dti/PD-343_visit_1_dti.nii.gz
for file in $(cat /home/neuro/codebank/GDPH/dti_list.txt); do
subj=$(basename $file .nii.gz)     ## HC-004_visit_2_dti
group_id=$(echo $subj | awk -F'_' '{print $1}' )
visit_id=$(echo $subj | awk -F'_' '{print $3}' )
echo "=============  subj is $subj ========="
tgt_dir=$(dirname $file)    ## /home/neuro/MRI_GDPH/nifti/HC-004/dti
echo "=============   tgt_dir is $tgt_dir =========="
# *dti.bval
cd $tgt_dir
mkdir -p temp

echo "===== Denoise and Degibbs ==== "

dwidenoise -nthreads 10 -force $file ./temp/"$subj"_temp_denoised.nii.gz
mrdegibbs  -nthreads 10 -force ./temp/"$subj"_temp_denoised.nii.gz ./temp/"$subj"_temp_denoised_degibbs.nii.gz

select_dwi_vols $file "$subj".bval ./temp/"$subj"_temp_b0s.nii.gz 0
bet ./temp/"$subj"_temp_b0s.nii.gz ./temp/"$subj"_temp_b0s_brain -f 0.3 -R -m

## eddy_
echo  "
====running eddy correction====
    "
eddy_openmp -v --imain=./temp/"$subj"_temp_denoised_degibbs.nii.gz --mask=./temp/"$subj"_temp_b0s_brain_mask.nii.gz  \
--acqp=/home/neuro/codebank/GDPH/philips_dti_acq_parameters.txt --index=/home/neuro/codebank/GDPH/philips_dti_index.txt \
--bvals="$subj".bval --bvecs="$subj".bvec --repol --out="$subj"-eddy_corrected_data

eddy_quad "$subj"-eddy_corrected_data -idx /home/neuro/codebank/GDPH/philips_dti_index.txt \
-par /home/neuro/codebank/GDPH/philips_dti_acq_parameters.txt -m ./temp/"$subj"_temp_b0s_brain_mask.nii.gz -b "$subj".bval

dtifit --verbose -k "$subj"-eddy_corrected_data.nii.gz -o "$subj"_dtifit \
-m ./temp/"$subj"_temp_b0s_brain_mask.nii.gz -r "$subj".bvec  -b "$subj".bval

echo "final brain mask"

select_dwi_vols "$subj"-eddy_corrected_data.nii.gz "$subj".bval "$subj"_temp_eddy_b0s.nii.gz 0
bet "$subj"_temp_eddy_b0s.nii.gz "$subj"_temp_eddy_b0s_brain -f 0.3 -R -m
fslmaths "$subj"_temp_eddy_b0s_brain_mask -fillh "$subj"-eddy_corrected_data-mask.nii.gz

echo "====== organize files  ====="

rm -r temp
rm *eddy_outlier_free_data.nii.gz
mkdir -p /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/eddy
mkdir -p /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/dtifit

 mv *.eddy_*   /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/eddy
 mv *_dtifit* /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/dtifit

mv *eddy_corrected_data.qc /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/
mv *eddy_corrected_data.nii.gz    /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/
mv *dti-eddy_corrected_data-mask.nii.gz /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/
mv *dti_temp_eddy_b0s.nii.gz /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/
mv *dti_temp_eddy_b0s_brain.nii.gz /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/
mv *dti_temp_eddy_b0s_brain_mask.nii.gz /home/neuro/MRI_GDPH/processed/"$group_id"/visit_"$visit_id"/dwi/

done

end=`date +%s`
runtime=$((end-start))
echo "runtime is $runtime"

6.2 DWI pre-processing with synthesized b0

#!/bin/bash

start=`date +%s`


cd /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti
nif_dir=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti

#/home/neuro/codebank/GDPH/ppmi_dwi_nif_list.txt
# "nif_dir"/3102/Visit-00/dwi/3102_Visit-00_dwi.nii.gz

# file=./14426/Visit-00/dwi/14426_Visit-00_dwi.nii.gz
for file in $(cat /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_dwi_nif_list.txt );
 do
  echo "file is $file"
subj=$(basename $file .nii.gz)     ## 3102_Visit-00_dwi
sub_id=$(echo $subj | awk -F'_' '{print $1}' )  ## 3102
visit_id=$(echo $subj | awk -F'_' '{print $2}' | awk -F'-' '{print $2}' )   ## 00

echo " ==== subj is is  $subj     ====
       ==== subj_id is $sub_id    ====
       ==== visit_id is $visit_id ==== "
#/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti/3104/Visit-00/dwi/3104_Visit-00_dwi.nii.gz
tgt_dir=$(dirname $file)   ## /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti/3104/Visit-00/dwi/
echo "=============   tgt_dir is $tgt_dir =========="

cd $tgt_dir
mkdir -p temp

echo "===== Denoise and Degibbs ==== "

dwidenoise -nthreads 10 -force $file ./temp/"$subj"_temp_denoised.nii.gz
mrdegibbs  -nthreads 10 -force ./temp/"$subj"_temp_denoised.nii.gz ./temp/"$subj"_temp_denoised_degibbs.nii.gz

select_dwi_vols $file "$subj".bval ./temp/"$subj"_temp_b0s.nii.gz 0
bet ./temp/"$subj"_temp_b0s.nii.gz ./temp/"$subj"_temp_b0s_brain -f 0.3 -R -m

end=`date +%s`
runtime=$((end-start))
echo "runtime is $runtime"

### synthesize b0

rm -r /home/neuro/INPUTS
rm -r /home/neuro/OUTPUTS
mkdir -p /home/neuro/INPUTS
mkdir -p /home/neuro/OUTPUTS

cp "$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/temp/"$sub_id"_Visit-"$visit_id"_dwi_temp_b0s.nii.gz  /home/neuro/INPUTS/b0.nii.gz
cp "$nif_dir"/"$sub_id"/Visit-"$visit_id"/t1/"$sub_id"_Visit-"$visit_id"_t1.nii.gz /home/neuro/INPUTS/T1.nii.gz
cp /home/neuro/codebank/GDPH/ppmi_dti_acq_parameters_synb0.txt /home/neuro/INPUTS/acqparams.txt

cd /home/neuro/
singularity run -e \
-B INPUTS/:/INPUTS \
-B OUTPUTS/:/OUTPUTS \
-B /usr/local/freesurfer/license.txt:/extra/freesurfer/license.txt /home/neuro/synb0-disco_v3.0.sif

end=`date +%s`
runtime=$((end-start))
echo "runtime is $runtime"

## eddy

echo  "
====running eddy correction====
   "
   ## /home/neuro/3102_Visit-00_dwi_eddy_corrected_data.nii.gz

eddy_openmp -v  --imain="$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/"$sub_id"_Visit-"$visit_id"_dwi.nii.gz \
--mask="$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/temp/"$sub_id"_Visit-"$visit_id"_dwi_temp_b0s_brain_mask.nii.gz \
--acqp=/home/neuro/codebank/GDPH/ppmi_dti_acq_parameters.txt \
--index=/home/neuro/codebank/GDPH/ppmi_dti_index.txt \
--bvecs="$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/"$sub_id"_Visit-"$visit_id"_dwi.bvec \
--bvals="$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/"$sub_id"_Visit-"$visit_id"_dwi.bval \
--topup=/home/neuro/OUTPUTS/topup  --repol --out="$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/"$subj"_eddy_corrected_data

echo "
 ===== QC for eddy =======
         "

cd $tgt_dir

# /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_HC_nifti/3102/Visit-00/dwi/3102_Visit-00_dwi_eddy_corrected_data.nii.gz

eddy_quad "$nif_dir"/"$sub_id"/Visit-"$visit_id"/dwi/"$subj"_eddy_corrected_data -idx /home/neuro/codebank/GDPH/ppmi_dti_index.txt \
-par /home/neuro/codebank/GDPH/ppmi_dti_acq_parameters.txt -m ./temp/"$subj"_temp_b0s_brain_mask.nii.gz -b "$subj".bval

dtifit --verbose -k "$subj"_eddy_corrected_data.nii.gz -o "$subj"_dtifit \
-m ./temp/"$subj"_temp_b0s_brain_mask.nii.gz -r "$subj".bvec  -b "$subj".bval

echo "==== final brain mask  ===="

select_dwi_vols "$subj"_eddy_corrected_data.nii.gz "$subj".bval "$subj"_eddy_b0s.nii.gz 0
bet "$subj"_eddy_b0s.nii.gz "$subj"_eddy_b0s_brain -f 0.3 -R -m
fslmaths "$subj"_eddy_b0s_brain_mask -fillh "$subj"_eddy_corrected_data-mask.nii.gz

echo "====== organize files  ====="

mkdir -p /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/eddy
mkdir -p /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/dtifit
mkdir -p /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/synb0_output
/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed
rm *eddy_outlier_free_data.nii.gz

 mv *.eddy_*   /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/eddy
 mv *_dtifit* /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/dtifit

mv *eddy_corrected_data.qc /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/
# mv *eddy_corrected_data.nii.gz    /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/
mv *dwi_eddy_corrected_data-mask.nii.gz /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/
mv *dwi_eddy_b0s.nii.gz /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/
mv *dwi_eddy_b0s_brain.nii.gz /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/
mv *dwi_eddy_b0s_brain_mask.nii.gz /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/


cp /home/neuro/OUTPUTS/b0_all.topup_log  /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/synb0_output/
cp /home/neuro/OUTPUTS/b0_all_topup.nii.gz /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/synb0_output/
cp /home/neuro/OUTPUTS/topup_fieldcoef.nii.gz /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/synb0_output/
cp /home/neuro/OUTPUTS/topup_movpar.txt /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/"$sub_id"/Visit-"$visit_id"/dwi/synb0_output/

rm -r ./temp/

done

end=`date +%s`
runtime=$((end-start))
echo "===========================
runtime is $runtime
## 3361s =55min
==========================="

6.3 PSMD calculation and extraction

#!/bin/bash

## dwi preprocessing

## Denoise and Degibbs

start=`date +%s`

PSMD=/home/neuro/codebank/psmd-main-20230126/psmd.sh
# /home/neuro/MRI_GDPH/processed/PD-277/visit_3/dwi/PD-277_visit_3_dti-eddy_corrected_data.nii.gz
# /home/neuro/MRI_GDPH/nifti/HC-064/dti/HC-064_visit_1_dti.bval
# /home/neuro/MRI_GDPH/nifti/HC-064/dti/HC-064_visit_1_dti.bvec
# /home/neuro/codebank/psmd-main-20230126/psmd.sh
# /home/neuro/codebank/psmd-main-20230126/skeleton_mask_2019.nii.gz
nif_dir=/home/neuro/MRI_GDPH/nifti
proc_dir=/home/neuro/MRI_GDPH/processed

# processed_dwi_list=$(find . -name '*eddy_corrected_data.nii*' | sort -n)
# echo "$processed_dwi_list" >> /home/neuro/codebank/GDPH/gdph_dwi_processed_list_20230126.txt

# dwi_file=/home/neuro/MRI_GDPH/processed/HC-004/visit_2/dwi/HC-004_visit_2_dti-eddy_corrected_data.nii.gz

for dwi_file in $(cat /home/neuro/codebank/GDPH/gdph_dwi_processed_list.txt);
do
subj=$(basename $dwi_file .nii.gz)     ## HC-004_visit_2_dti-eddy_corrected_data
echo  "$subj"
group_id=$(echo $subj | awk -F'_' '{print $1}' )  ## HC-004
visit_id=$(echo $subj | awk -F'_' '{print $3}' )   ## 2
echo "=============  subj is $subj
  group_id is $group_id
  visit_id is $visit_id  ========="

bval_file=/home/neuro/MRI_GDPH/nifti/$group_id/dti/"$group_id"_visit_"$visit_id"_dti.bval
bvec_file=/home/neuro/MRI_GDPH/nifti/$group_id/dti/"$group_id"_visit_"$visit_id"_dti.bvec

mkdir -p /home/neuro/MRI_GDPH/psmd_metrics_gdph/temp_"$subj"
cd /home/neuro/MRI_GDPH/psmd_metrics_gdph/temp_"$subj"

echo "====Runing PSMD.....====="
$PSMD -d $dwi_file -b $bval_file -r $bvec_file -s /home/neuro/codebank/psmd-main-20230126/skeleton_mask_2019.nii.gz  \
> /home/neuro/MRI_GDPH/psmd_metrics_gdph/PSMD_"$subj".txt

done
# rm *.nii.gz # remove dwi images
# cat /project/3015083.01/mengfei/PSMD_2011/PSMD_001/PSMD_001.txt | awk 'FNR==2 {print}'
# /home/neuro/MRI_GDPH/psmd_metrics_gdph

# ## print 11th row, 3rd column,-->psmd value

# ID=/home/neuro/MRI_GDPH/psmd_metrics_gdph/PSMD_HC-004_visit_2_dti-eddy_corrected_data.txt
for ID in $(cat /home/neuro/codebank/GDPH/gdph_psmd_list.txt);
do
ID_easy=$(echo $ID | awk -F'PSMD_' '{print $2}' | awk -F'_dti' '{print $1}' )
echo "ID easy is $ID_easy"
psmd_value=$(cat $ID | awk 'FNR==11 {print}' | awk -F" " '{print $3}' )
echo "psmd is $psmd_value"
cd /home/neuro/codebank/GDPH
echo $ID_easy $psmd_value >> gdph_psmd_values_20230129-final.txt
#gdph_psmd_values_20230129.txt
done

#
end=`date +%s`
runtime=$((end-start))
echo "runtime is $runtime"

6.4 Extract 3D volumes from 4D volumes

# Environment
#module unload fsl

#module load fsl/6.0.1

niftifolder=/project/3015083.01/nifti

cd "${niftifolder}"  # this is home folder for nifti folders

## for subject in $(ls | head -3);
# for subject in $(ls);

subject="$1"
echo "subject ID: $subject"
qsmfolder="${niftifolder}"/"${subject}"/ses-mri01/qsm

echo " ok 1" 

if [ -d "${qsmfolder}" ]; then 
    cd "${qsmfolder}" 
    echo "current folder: $PWD"
    if [ -f "${subject}"_ses-mri01_qsm_mag9.nii.gz ]
    then
        echo "let's skip"
    else        
        echo "Moi, let's do the splitting..."
# extract each magnitude volume from 4D data
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag1 --vols=0
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag2 --vols=1
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag3 --vols=2

    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag4 --vols=3
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag5 --vols=4
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag6 --vols=5

    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag7 --vols=6
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag8 --vols=7
    fslselectvols -i "${subject}"_ses-mri01_qsm.nii.gz -o "${subject}"_ses-mri01_qsm_mag9 --vols=8
    echo "split of mag data was done"

# extract each phase volume from 4D data
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph1 --vols=0
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph2 --vols=1
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph3 --vols=2

    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph4 --vols=3
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph5 --vols=4
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph6 --vols=5

    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph7 --vols=6
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph8 --vols=7
    fslselectvols -i "${subject}"_ses-mri01_qsm_ph.nii.gz -o "${subject}"_ses-mri01_qsm_ph9 --vols=8
    echo "split of phase data was done"
    fi
fi

6.5 Tractseg pipeline

ppmi_dir=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi
##file=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti/3390/Visit-04/dwi/3390_Visit-04_dwi_eddy_corrected_data.nii.gz

for file in  $(cat /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_tract_seg.txt);
do
  start=`date +%s`
  echo "==========================
  dwi file is $file"
group_name=$(echo $file | awk -F'/' '{print $7}' ) ## ppmi_HC_prodromal_nifti or ppmi_PD_nifti
group_name_processed=$(echo $group_name | sed 's/nifti/processed/')
subj_id=$(echo $file | awk -F'/' '{print $8}' ) ## 3390
visit_id=$(echo $file | awk -F'/' '{print $9}' )  ## Visit-04

echo "=====================
group_name is $group_name
subj_id is $subj_id
visit_id is  $visit_id
group_name_processed is  $group_name_processed"


echo "=============  subj is $subj_id ========="

# mkdir -p "$tgt_dir"/tractseg
echo "=============   tgt_dir is $tgt_dir =========="

tgt_dir=$(dirname $file)    ## /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti/3390/Visit-04/dwi

mkdir -p "$tgt_dir"/tractseg
cd $tgt_dir/tractseg
echo "========  PWD is $PWD  ========"

cp "$tgt_dir"/"$subj_id"_"$visit_id"_dwi_eddy_corrected_data.nii.gz dwi.nii.gz
cp "$tgt_dir"/"$subj_id"_"$visit_id"_dwi.bval dwi.bval
cp "$tgt_dir"/"$subj_id"_"$visit_id"_dwi.bvec dwi.bvec

echo "============= TractSeg is running ========="

TractSeg -i dwi.nii.gz -o . --bvals dwi.bval --bvecs dwi.bvec --raw_diffusion_input

rm dwi*
mv "$tgt_dir"/tractseg  /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/

# /hsome/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_processed/14426/Visit-00/dwi

cd /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations

echo "=====merge L+R masks======="
fslmaths AF_left.nii.gz -add AF_right.nii.gz   all_AF.nii.gz
fslmaths ATR_left.nii.gz -add ATR_right.nii.gz   all_ATR.nii.gz
fslmaths CG_left.nii.gz -add CG_right.nii.gz   all_CG.nii.gz
fslmaths CST_left.nii.gz -add CST_right.nii.gz   all_CST.nii.gz
fslmaths FPT_left.nii.gz -add FPT_right.nii.gz   all_FPT.nii.gz
fslmaths FX_left.nii.gz -add FX_right.nii.gz   all_FX.nii.gz
fslmaths ICP_left.nii.gz -add ICP_right.nii.gz   all_ICP.nii.gz
fslmaths IFO_left.nii.gz -add IFO_right.nii.gz   all_IFO.nii.gz
fslmaths ILF_left.nii.gz -add ILF_right.nii.gz   all_ILF.nii.gz
fslmaths MLF_left.nii.gz -add MLF_right.nii.gz   all_MLF.nii.gz
fslmaths OR_left.nii.gz -add OR_right.nii.gz   all_OR.nii.gz
fslmaths POPT_left.nii.gz -add POPT_right.nii.gz   all_POPT.nii.gz
fslmaths SCP_left.nii.gz -add SCP_right.nii.gz   all_SCP.nii.gz
fslmaths SLF_III_left.nii.gz -add SLF_III_right.nii.gz   all_SLF_III.nii.gz
fslmaths SLF_II_left.nii.gz -add SLF_II_right.nii.gz   all_SLF_II.nii.gz
fslmaths SLF_I_left.nii.gz -add SLF_I_right.nii.gz   all_SLF_I.nii.gz
fslmaths ST_FO_left.nii.gz -add ST_FO_right.nii.gz   all_ST_FO.nii.gz
fslmaths ST_OCC_left.nii.gz -add ST_OCC_right.nii.gz   all_ST_OCC.nii.gz
fslmaths ST_PAR_left.nii.gz -add ST_PAR_right.nii.gz   all_ST_PAR.nii.gz
fslmaths ST_POSTC_left.nii.gz -add ST_POSTC_right.nii.gz   all_ST_POSTC.nii.gz
fslmaths ST_PREC_left.nii.gz -add ST_PREC_right.nii.gz   all_ST_PREC.nii.gz
fslmaths ST_PREF_left.nii.gz -add ST_PREF_right.nii.gz   all_ST_PREF.nii.gz
fslmaths ST_PREM_left.nii.gz -add ST_PREM_right.nii.gz   all_ST_PREM.nii.gz
fslmaths STR_left.nii.gz -add STR_right.nii.gz   all_STR.nii.gz
fslmaths T_OCC_left.nii.gz -add T_OCC_right.nii.gz   all_T_OCC.nii.gz
fslmaths T_PAR_left.nii.gz -add T_PAR_right.nii.gz   all_T_PAR.nii.gz
fslmaths T_POSTC_left.nii.gz -add T_POSTC_right.nii.gz   all_POSTC_PAR.nii.gz
fslmaths T_PREC_left.nii.gz -add T_PREC_right.nii.gz   all_PREC_PAR.nii.gz
fslmaths T_PREF_left.nii.gz -add T_PREF_right.nii.gz   all_PREF_PAR.nii.gz
fslmaths T_PREM_left.nii.gz -add T_PREM_right.nii.gz   all_PREM_PAR.nii.gz
fslmaths UF_left.nii.gz -add UF_right.nii.gz   all_UF.nii.gz

mkdir -p /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations_merged
mv  /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations/all*nii.gz \
/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations_merged

end=`date +%s`
runtime=$((end-start))
echo "runtime is $runtime"
done

6.6 Tractseg FA_MD metrics extraction

#!/bin/bash

# find -name  *dwi_eddy_corrected_data.nii.gz | sort -n > ./ppmi_tract_seg_1.txt
# /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti/3106/Visit-00/dwi/3106_Visit-00_dwi_eddy_corrected_data.nii.gz
#
# ./ppmi_HC_prodromal_nifti/3390/Visit-04/dwi/3390_Visit-04_dwi_eddy_corrected_data.nii.gz
#

ppmi_dir=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi
##file=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_HC_prodromal_nifti/3390/Visit-04/dwi/3390_Visit-04_dwi_eddy_corrected_data.nii.gz

for file in  $(cat /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_tract_seg.txt );
do
  start=`date +%s`
  echo "==========================
  dwi file is $file"
group_name=$(echo $file | awk -F'/' '{print $7}' ) ## ppmi_HC_prodromal_nifti or ppmi_PD_nifti
group_name_processed=$(echo $group_name | sed 's/nifti/processed/')   ## ppmi_HC_prodromal_processed or ppmi_PD_processed
subj_id=$(echo $file | awk -F'/' '{print $8}' ) ## 3390
visit_id=$(echo $file | awk -F'/' '{print $9}' )  ## Visit-04

echo "=====================
group_name is $group_name
subj_id is $subj_id
visit_id is  $visit_id
group_name_processed is  $group_name_processed"

      md_file=/"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/dtifit/"$subj_id"_"$visit_id"_dwi_dtifit_MD.nii.gz
      fa_file=/"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/dtifit/"$subj_id"_"$visit_id"_dwi_dtifit_FA.nii.gz
      roi_dir=/"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations_merged

      cp /"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations/CA.nii.gz /"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations_merged/all_CA.nii.gz
      cp /"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations/CC.nii.gz /"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations_merged/all_CC.nii.gz
      cp /"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations/MCP.nii.gz /"$ppmi_dir"/"$group_name_processed"/"$subj_id"/"$visit_id"/dwi/tractseg/bundle_segmentations_merged/all_MCP.nii.gz

      for roi_file in ${roi_dir}/all*.nii.gz; do
      roi_name=$(basename ${roi_file} .nii.gz )
      md_mean=$(fslstats ${md_file} -k ${roi_file} -M )
      fa_mean=$(fslstats ${fa_file} -k ${roi_file} -M )
      echo "
      =================
      roi_name: MD_mean FA_mean == ${roi_name}: ${md_mean} ${fa_mean}
      =================
      "
      echo "ID_visit_group tract_name: MD_mean FA_mean == "$subj_id"_"$visit_id"_"${roi_name}"_"$group_name_processed": ${md_mean} ${fa_mean}" >> /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_tractseg_FA_MD_metrics_20230427.txt
      done

end=`date +%s`
runtime=$((end-start))
echo "runtime is $runtime"
done

##============================================

6.7 Freewater pipeline

# Path to the file containing the list of text file paths
file_list="/home/neuro/caixin/dwi_corrected_list_0520.txt"

# /home/neuro/caixin/processed/HC-001/ses-1/dwi/HC-001_eddy_corrected_data.nii.gz

# Read each line from the file list
for dwi in $(cat "$file_list" )
do
    echo "=========================Processing file: $dwi"
    start=`date +%s`
    #dwi # /home/neuro/caixin/processed/HC-001/ses-1/dwi/HC-001_eddy_corrected_data.nii.gz

    group_id=$(echo "$dwi"  | awk -F'/' '{print $6}' )  # HC-001
    echo "====== $group_id ======="

  

in_dti=/home/neuro/caixin/processed/"$group_id"/ses-1/dwi/"$group_id"_eddy_corrected_data.nii.gz

in_dti_mask=/home/neuro/caixin/processed/"$group_id"/ses-1/dwi/"$group_id"_eddy_b0s_brain_mask.nii.gz
in_bvec=/home/neuro/caixin/nifti/"$group_id"/ses-1/dti/"$group_id"_dti.bvec

in_bval=/home/neuro/caixin/processed/file2.bval  # change to b0 and b1000 is necessary to get FW

  mkdir -p  /home/neuro/caixin/processed/"$group_id"/ses-1/freewater
 
  cd  /home/neuro/caixin/processed/"$group_id"/ses-1/freewater  

/home/neuro/codebank/scripts_FW_CONSORTIUM_free_water/MAIN_script_FW.sh ${in_dti} ${in_dti_mask} ${in_bval} ${in_bvec} .

done

6.8 Free water in T1 longitudinal space

export SUBJECTS_DIR=/media/neuro/16t/ppmi_t1_long_all/t1_longitudinal_freesurfer_template_all

# /media/neuro/16t/ppmi_t1_long_all/3102_v06.long.3102_template
# /media/neuro/16t/ppmi_t1_long_all/3105_BL.long.3105_template

for subject in $(cat /home/neuro/codebank/GDPH/ppmi/t1_longitudinal_freesurfer_template_all.txt);
do 
echo "=======sub is  $subject ====="

# 3102_BL.long.3102_template
# 3102_v06.long.3102_template

id=$(echo $subject | awk -F'_' '{print $1}' )  ## 3102
visit_id=$(echo $subject | awk -F'.' '{print $1}' | awk -F'_' '{print $2}' )  #3102_BL --> BL


if [ "$visit_id" == "BL" ]; then
    visit_id="00"
elif [ "$visit_id" == "v04" ]; then
    visit_id="04"
elif [ "$visit_id" == "v06" ]; then
    visit_id="06"
elif [ "$visit_id" == "v10" ]; then
    visit_id="10"
fi

echo "==========subject, id ,and Converted visit_id:  $subject 
 $id 
  $visit_id=========="

FA_image=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/"$id"/Visit-"$visit_id"/dwi/dtifit/"$id"_Visit-"$visit_id"_dwi_dtifit_FA.nii.gz 
MD_image=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/"$id"/Visit-"$visit_id"/dwi/dtifit/"$id"_Visit-"$visit_id"_dwi_dtifit_MD.nii.gz
FW_image=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/"$id"/Visit-"$visit_id"/dwi/freewater/wls_dti_FW.nii.gz
MDt_image=/home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/"$id"/Visit-"$visit_id"/dwi/freewater/fwc_wls_dti_MD.nii.gz

# /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/3105/Visit-06/dwi/freewater/fwc_wls_dti_FA.nii.gz 
# /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/3105/Visit-06/dwi/freewater/fwc_wls_dti_MD.nii.gz 
# /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/3105/Visit-06/dwi/freewater/wls_dti_FW.nii.gz

    FA_in_fs_space=${FA_image%.nii.gz}_fs.nii.gz
    MD_in_fs_space=${MD_image%.nii.gz}_fs.nii.gz
    FW_in_fs_space=${FW_image%.nii.gz}_fs.nii.gz
    MDt_in_fs_space=${MDt_image%.nii.gz}_fs.nii.gz

    mkdir -p /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/"$id"/Visit-"$visit_id"/dwi/freesurfer_longitudinal
    cd /home/neuro/MRI_GDPH/ppmi/20220919_ppmi/ppmi_PD_processed/"$id"/Visit-"$visit_id"/dwi/freesurfer_longitudinal

    bbregister --s $subject --mov $FA_image --init-fsl --reg register.dat --o $FA_in_fs_space --dti   # to get register.dat
    mri_vol2vol --mov $FW_image --targ $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --o $FW_in_fs_space --reg register.dat --interp trilin
    mri_vol2vol --mov $MDt_image --targ $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --o $MDt_in_fs_space --reg register.dat --interp trilin
    
 echo "     
   ========= current dir is $PWD =======
 "

    for hemi in lh rh; do
    mri_vol2surf --mov $FW_in_fs_space --hemi $hemi --surf white --projfrac 0.5 \
    --o ${subject}_${hemi}_FW.mgh --reg register.dat --trgsubject $subject --surf-fwhm 10 --cortex
    done


   for hemi in lh rh; do
    mri_vol2surf --mov $MDt_in_fs_space --hemi $hemi --surf white --projfrac 0.5 \
    --o ${subject}_${hemi}_MDt.mgh --reg register.dat --trgsubject $subject --surf-fwhm 10 --cortex
    done

    for hemi in lh rh; do
    mri_segstats --annot ${subject} ${hemi} aparc --in ${subject}_${hemi}_FW.mgh --sum ${subject}_${hemi}.fw.sum
    done

    for hemi in lh rh; do
    mri_segstats --annot ${subject} ${hemi} aparc --in ${subject}_${hemi}_MDt.mgh --sum ${subject}_${hemi}.MDt.sum
    done


done