Update stats plots, add longitudinal sample size calculation

.gitignore

@@ -10,3 +10,5 @@ qc/
 *.jpg
 *.png
 .gitignore
+*.egg-info/
+dist/

README.md

@@ -69,9 +69,10 @@ Per-subject stat: Panda dataframe `df_sub`:
 - intra-subject MEAN: MEAN[CSA(sI, rX, :)] --> MEAN_intra(sI, rX): `df_sub['mean']`
 - intra-subject STD: STD[CSA(sI, rX, :)] --> STD_intra(sI, rX): `df_sub['std']`
 - intra-subject COV: STD[CSA(sI, rX, :)] / MEAN[CSA(sI, rX, :)] --> COV_intra(sI, rX): `df_sub['cov']`
-- rescale_estimated_subject MEAN: MEAN[CSA(sI, rX, :) / CSA(sI, 1, :)] --> MEAN_rescale_estimated_subject(sI, rX): `df_sub['rescale_estimated']`
+- rescale_estimated_subject MEAN: MEAN[ CSA(sI, rX, :) ] / MEAN[ CSA(sI, 1, :) ] --> MEAN_rescale_estimated_subject(sI, rX): `df_sub['rescale_estimated']`
 - intra-subject error MEAN: MEAN[CSA(sI, rX, :)] - (rX^2 * MEAN[CSA(sI, 1, :)]) --> MEAN_error_intra(sI, rX): `df_sub['error']`
 - intra-subject error in percentage MEAN: [MEAN[CSA(sI, rX, :)] - (rX^2 * MEAN[CSA(sI, 1, :)])] / MEAN[CSA(sI, rX, :)] --> MEAN_error_intra_perc(sI, rX): `df_sub['perc_error']`
+- intra-subject difference MEAN: [MEAN[CSA(sI, 1, :)] - (rX^2 * MEAN[CSA(sI, rX, :)])] --> MEAN_diff(sI, rX): `df_sub['diff']`
 
 Across-subject stats: Panda dataframe `df_rescale`
 - intra-subject STD: MEAN[STD_intra(:, rX)] --> STD_intra(rX): `df_rescale['std_intra']`
@@ -81,9 +82,11 @@ Across-subject stats: Panda dataframe `df_rescale`
 - rescale_estimated (across subjects) STD: STD[MEAN_rescale_estimated_subject(:, rX)] --> STD_rescale_estimated(rX): `df_rescale['std_rescale_estimated']`
 - error in percentage (across subjects) MEAN: MEAN[MEAN_error_intra(:, rX)]
 - error in percentage (across subjects) STD: STD[MEAN_error_intra(:, rX)]
+- intra-subject difference STD: STD[MEAN_diff(:, rX)]
 
 Power analysis:
-- sample size: [(z(uncertainty) + z(power))^2 * (2 * STD[MEAN(:, rX)]^2)] / [MEAN[CSA(sI, 1, :)] - MEAN[CSA(sI, rX, :)]] 
+- sample size for a between subject study: [(z(uncertainty) + z(power))^2 * (2 * STD[MEAN(:, rX)]^2)] / [MEAN[CSA(sI, 1, :)] - MEAN[CSA(sI, rX, :)]] 
+- sample size for a within subject study: [(z(uncertainty) + z(power))^2 * (STD[MEAN_diff(:, rX)]^2)] / [MEAN[CSA(sI, 1, :)] - MEAN[CSA(sI, rX, :)]] 
 
 Plot results:
 - STD_intersub

affine_transfo.py

@@ -120,7 +120,11 @@ def random_values(df, subject_name, config_param):
         'shift_PA': shift_PA,
         'shift_IS': shift_IS
     }
-    df = df.append(transfo_dict, ignore_index=True)
+    #df = df.append(transfo_dict, ignore_index=True)
+    print(df)
+    print(transfo_dict)
+    df_new_row = pd.DataFrame(transfo_dict, index=[0])
+    df = pd.concat([df, df_new_row], ignore_index=True)#.reset_index()
     return df, angle_IS, angle_PA, angle_LR, shift_LR, shift_PA, shift_IS
 
 

config_sct_run_batch.yml

@@ -1,10 +1,9 @@
 # config file for sct_run_batch
-path_data: /scratch/pabaua/data-multi-subject-p
-path_output: /scratch/pabaua/results_csa_t2
-script: /home/pabaua/csa-atrophy/process_data.sh
-script_args: /home/pabaua/csa-atrophy/config_script.yml
-jobs: -1
+path_data: ~/datasets/data-multi-subject
+path_output: ~/process_results/results_csa-atrophy_t2-all-2024-08-13
+script: ~/code/csa-atrophy/process_data.sh
+script_args: /home/GRAMES.POLYMTL.CA/sebeda/code/csa-atrophy/config_script.yml
+jobs: 8
 batch_log: sct_run_batch_log.txt
 continue_on_error: 1
 subject_prefix: sub-
-zip: true

config_sct_run_batch_cc.yml

@@ -0,0 +1,10 @@
+# config file for sct_run_batch
+path_data: /home/sabeda/scratch/data-multi-subject
+path_output: /home/sabeda/scratch/results_csa-atrophy_t2-all
+script: /home/sabeda/csa-atrophy/process_data.sh
+script_args: /home/sabeda/csa-atrophy/config_script.yml
+jobs: -1
+batch_log: sct_run_batch_log.txt
+continue_on_error: 1
+subject_prefix: sub-
+zip: true

job_template.sh

@@ -3,7 +3,7 @@
 #SBATCH --nodes=1
 #SBATCH --cpus-per-task=32    # number of OpenMP processes
 #SBATCH --mem=128G
-#SBATCH --mail-user=paul.bautin@polymtl.ca
+#SBATCH --mail-user=sandrine.bedard@polymtl.ca
 #SBATCH --mail-type=ALL
 export OMP_NUM_THREADS=1
 export MKL_NUM_THREADS=1

process_data.sh

@@ -49,58 +49,63 @@ segment_and_label_if_does_not_exist(){
   local file="$1"
   local contrast="$2"
   local contrast_str="$3"
-  segment_if_does_not_exist $file ${contrast} "-qc ${PATH_QC} -qc-subject ${SUBJECT}"
+  FILESEG="${file}_seg"
+  sct_deepseg -i ${file}.nii.gz -task seg_sc_contrast_agnostic #-qc $PATH_QC -qc-subject ${SUBJECT}
+  #sct_qc -i ${file}.nii.gz -s ${FILESEG}.nii.gz -p sct_deepseg_sc -qc $PATH_QC -qc-subject ${SUBJECT}
+  #segment_if_does_not_exist $file ${contrast} "-qc ${PATH_QC} -qc-subject ${SUBJECT}"
   local file_seg=${FILESEG}
   # Update global variable with segmentation file name
-  FILELABEL="${file}_labels-disc"
-  FILELABELMANUAL="${path_derivatives}/${SUBJECT}_${contrast_str}_labels-disc-manual"
+  FILELABEL="${file}_label-discs_dlabel"  # TODO change for updated names
+  FILELABELMANUAL="${path_derivatives}/${SUBJECT}_${contrast_str}_label-discs_dlabel"
   if [ -e "${FILELABELMANUAL}.nii.gz" ]; then
     echo "manual labeled file was found: ${FILELABELMANUAL}"
+    rsync -avzh ${FILELABELMANUAL}.nii.gz ${FILELABEL}.nii.gz
     # reorienting and resampling image
-    sct_image -i ${FILELABELMANUAL}.nii.gz -setorient RPI -o "${FILELABELMANUAL}_RPI.nii.gz"
-    sct_maths -i ${FILELABELMANUAL}_RPI.nii.gz -dilate 2 -o ${FILELABELMANUAL}_RPI_dil.nii.gz
-    sct_resample -i ${FILELABELMANUAL}_RPI_dil.nii.gz -mm $interp -x nn -o ${FILELABELMANUAL}_RPI_dil_r.nii.gz
-    rsync -avzh "${FILELABELMANUAL}_RPI_dil_r.nii.gz" ${FILELABEL}.nii.gz
+    #sct_image -i ${FILELABELMANUAL}.nii.gz -setorient RPI -o "${FILELABELMANUAL}_RPI.nii.gz"
+    #sct_maths -i ${FILELABELMANUAL}_RPI.nii.gz -dilate 2 -o ${FILELABELMANUAL}_RPI_dil.nii.gz
+    #sct_resample -i ${FILELABELMANUAL}_RPI_dil.nii.gz -mm $interp -x nn -o ${FILELABELMANUAL}_RPI_dil_r.nii.gz
+    #rsync -avzh "${FILELABELMANUAL}_RPI_dil_r.nii.gz" ${FILELABEL}.nii.gz
     # Generate labeled segmentation
-    sct_label_vertebrae -i ${file}.nii.gz -s ${file_seg}.nii.gz -c ${contrast} -discfile "${FILELABELMANUAL}_RPI_dil_r.nii.gz" -qc ${PATH_QC} -qc-subject ${SUBJECT}
+    sct_image -i ${file_seg}.nii.gz -set-sform-to-qform
+    sct_image -i ${file}.nii.gz -set-sform-to-qform
+    sct_image -i "${FILELABEL}.nii.gz" -set-sform-to-qform
+    sct_label_vertebrae -i ${file}.nii.gz -s ${file_seg}.nii.gz -c ${contrast} -discfile "${FILELABEL}.nii.gz" #-qc ${PATH_QC} -qc-subject ${SUBJECT}
   else
     # Generate labeled segmentation
-    sct_label_vertebrae -i ${file}.nii.gz -s ${file_seg}.nii.gz -c ${contrast} -qc ${PATH_QC} -qc-subject ${SUBJECT}
+    sct_label_vertebrae -i ${file}.nii.gz -s ${file_seg}.nii.gz -c ${contrast} #-qc ${PATH_QC} -qc-subject ${SUBJECT}
   fi
   # Create labels in the Spinal Cord
-  sct_label_utils -i ${file_seg}_labeled.nii.gz -vert-body 0 -o ${FILELABEL}.nii.gz
+  #sct_label_utils -i ${file_seg}_labeled.nii.gz -vert-body 0 -o ${FILELABEL}.nii.gz
   FILE_SEG_LABEL=${file_seg}_labeled
 }
 
 # Check if manual segmentation already exists. If it does, copy it locally. If
 # it does not, perform seg.
-segment_if_does_not_exist(){
-  local file="$1"
-  local contrast="$2"
-  local qc=$3
-  # Update global variable with segmentation file name
-  FILESEG="${file}_seg"
-  FILESEGMANUAL="${path_derivatives}/${FILESEG}-manual"
-  if [ -e "${FILESEGMANUAL}.nii.gz" ]; then
-    echo "Found! Using manual segmentation."
-    sct_resample -i ${FILESEGMANUAL}.nii.gz -mm $interp -x nn -o ${FILESEGMANUAL}_r.nii.gz
-    rsync -avzh ${FILESEGMANUAL}_r.nii.gz ${FILESEG}.nii.gz
-    sct_qc -i ${file}.nii.gz -s ${FILESEG}.nii.gz -p sct_deepseg_sc $qc
-  else
-    # Segment spinal cord
-    sct_deepseg_sc -i ${file}.nii.gz -c ${contrast} $qc
-  fi
-}
+# segment_if_does_not_exist(){
+#   local file="$1"
+#   local contrast="$2"
+#   local qc=$3
+#   # Update global variable with segmentation file name
+#   FILESEG="${file}_seg"
+#   FILESEGMANUAL="${path_derivatives}/${FILESEG}-manual"
+#   if [ -e "${FILESEGMANUAL}.nii.gz" ]; then
+#     echo "Found! Using manual segmentation."
+#     sct_resample -i ${FILESEGMANUAL}.nii.gz -mm $interp -x nn -o ${FILESEGMANUAL}_r.nii.gz
+#     rsync -avzh ${FILESEGMANUAL}_r.nii.gz ${FILESEG}.nii.gz
+#     sct_qc -i ${file}.nii.gz -s ${FILESEG}.nii.gz -p sct_deepseg_sc $qc
+#   else
+#     # Segment spinal cord
+#     sct_deepseg_sc -i ${file}.nii.gz -c ${contrast} $qc
+#   fi
+# }
 
 
 # SCRIPT STARTS HERE
 # ==============================================================================
 # Display useful info for the log, such as SCT version, RAM and CPU cores available
 sct_check_dependencies -short
 # copy derivatives directory containing manual corrections to PATH_DATA_PROCESSED
-mkdir -p "${PATH_DATA_PROCESSED}/derivatives/labels/${SUBJECT}/anat/"
-cp -r "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat" "${PATH_DATA_PROCESSED}/derivatives/labels/${SUBJECT}"
-path_derivatives="${PATH_DATA_PROCESSED}/derivatives/labels/${SUBJECT}/anat"
+path_derivatives="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat"
 # Go to results folder, where most of the outputs will be located
 cd $PATH_DATA_PROCESSED
 # Copy source images
@@ -124,13 +129,13 @@ fi
 file=${SUBJECT}_${contrast_str}
 # Reorient image to RPI
 sct_image -i ${file}.nii.gz -setorient RPI -o ${file}_RPI.nii.gz
-file=${file}_RPI
 # Resample image isotropically
-sct_resample -i ${file}.nii.gz -mm $interp -o ${file}_r.nii.gz
-file=${file}_r
-end=`date +%s`
-runtime=$((end-start))
-echo "+++++++++++ TIME: Duration of reorienting and resampling:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
+sct_resample -i ${file}_RPI.nii.gz -mm $interp -o ${file}_RPI_r.nii.gz
+mv "${file}_RPI_r.nii.gz" "${SUBJECT}_space-other_${contrast_str}.nii.gz"
+file=${SUBJECT}_space-other_${contrast_str}
+# end=`date +%s`
+# runtime=$((end-start))
+# echo "+++++++++++ TIME: Duration of reorienting and resampling:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
 
 # Label spinal cord (only if it does not exist) in dir anat
 start=`date +%s`
@@ -198,9 +203,13 @@ for r_coef in ${R_COEFS[@]}; do
     file_label_r_t=${file_label_r_t}_crop
 
 
-    # Segment spinal cord (only if it does not exist)
+    # Segment spinal cord
     start=`date +%s`
-    sct_deepseg_sc -i ${file_r_t}.nii.gz -c ${contrast}
+    #sct_deepseg_sc -i ${file_r_t}.nii.gz -c ${contrast}
+    sct_deepseg -i ${file_r_t}.nii.gz -task seg_sc_contrast_agnostic # -qc $PATH_QC -qc-subject ${SUBJECT}
+    # TODO: soft csa?
+    #sct_qc -i ${file_r_t}.nii.gz -s ${file_r_t}_seg.nii.gz -p sct_deepseg_sc -qc $PATH_QC -qc-subject ${SUBJECT}
+
     end=`date +%s`
     runtime=$((end-start))
     echo "+++++++++++ TIME: Duration of segmentation t${i_transfo}:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"

requirements.txt → requirements_cc.txt

@@ -5,9 +5,9 @@ matplotlib
 scikit-image
 nibabel
 ruamel.yaml
+argparse
+setuptools
+PyYAML
+coloredlogs
+networkx
 yaml-1.3
-argparse~=1.4.0
-setuptools~=49.6.0
-PyYAML~=5.3.1
-coloredlogs~=14.0
-

requirements_manual_pip.txt

@@ -0,0 +1,2 @@
+networkx==3.2.1
+yaml-1.3

setup.py

@@ -29,6 +29,7 @@
     ],
     keywords='',
     install_requires=install_reqs,
+    packages = [],
     entry_points={
         'console_scripts': [
             'affine_transfo=affine_transfo:main',