localization-VBM.sh

#!/bin/sh

clear


     gm4D="/Users/stowler/ZZ.VBM.2014/inputs/GM_mod_merg_s2.nii.gz"
  blockLH="/Users/stowler/ZZ.VBM.2014/inputs/MNI152_T1_1mm_blockMask_LH.nii.gz"
  blockRH="/Users/stowler/ZZ.VBM.2014/inputs/MNI152_T1_1mm_blockMask_RH.nii.gz"
 idMatrix="/Users/stowler/ZZ.VBM.2014/inputs/id.mat"
parentDir="/Users/stowler/ZZ.VBM.2014/masked.GM_thresh.z165"


# list of models represented as their existing subdirectory names:
modelList="zzReq01_age+fit+exec+age.fit+exec.fit+.feat zzReq02_age+pa+exec+age.pa+exec.pa+.feat zzReq03_age+fitness+pa+age.fit+age.pa+.feat zzReq04_age+fit+age.fit+.feat zzReq05_age+pa+age.pa.feat"
echo ""
echo "Each VBM model has a directory of results in which there are masks that we are"
echo "going to localize. The directories are named for the models:"
echo ""
for model in ${modelList}; do
	ls -ald ${parentDir}/${model}
done


# loop start: perform cluster localization for clusters within each model:
for model in ${modelList}; do

	echo ""
	echo "############################################################"
	echo "# MODEL: ${parentDir}/${model}"
	echo "############################################################"
	echo ""

	cd ${parentDir}/${model}

	# loop start: each file containing cluster masks will be processed in its own iteration of this loop:
	for clusterMask in cluster_mask_zstat?.nii.gz; do
		echo ""
		echo "############################################################"
		echo "# MODEL: ${parentDir}/${model}"
		echo "#     CLUSTER MASK: ${clusterMask}"
		echo "############################################################"
		echo ""
		echo "Cluster mask:"
		ls -l ${clusterMask}
		fslinfo ${clusterMask}
		#rm -f ${parentDir}/meanGM-${model}++++${clusterMask}.1D
		#3dROIstats -mask ${clusterMask} ${gm4D} > ${parentDir}/meanGM-${model}++++${clusterMask}.1D

	
		# (TBD below: generalize this by wrapping in a test: does input cluster match 1mm MNI152 geometry already?)
		#
		echo ""
		echo "Converting from 2mm MNI152 resolution (odd number of sagittal"
		echo "slices) to 1mm MNI 152 resolution (even number of sagittal"
		echo "slices) so that the volume can be split into left and right halves"
		echo "of equal numbers of saggital slices:"
		rm -f ${clusterMask}-1mm.nii.gz
		flirt \
		-interp nearestneighbour \
		-datatype char \
		-in ${clusterMask} \
		-ref $FSLDIR/data/standard/MNI152_T1_1mm.nii.gz \
		-init ${idMatrix} \
		-applyisoxfm 1 \
		-out ${clusterMask}-1mm 
		echo "...done:"
		ls -altr ${clusterMask}-1mm.nii.gz
		#fslinfo ${clusterMask}-1mm.nii.gz
		#3dinfo ${clusterMask}-1mm.nii.gz
		
		# provide some evidence that upsampling didn't destroy the input mask:
		# (TBD maybe add more evidence, like histogram/mean/sd/etc)
		echo ""
		echo "These are the unique non-zero intensities in the upsampled mask per 3dmaskdump:"
		maskValues=`3dmaskdump -noijk -nozero -quiet ${clusterMask}-1mm.nii.gz | sort | uniq`
		echo "${maskValues}"
		echo ""
		echo "...they should be identical to the unique non-zero intensities in original mask per 3dmaskdump:"
		maskValuesLowRes=`3dmaskdump -noijk -nozero -quiet ${clusterMask} | sort | uniq`
		echo "${maskValuesLowRes}"
		echo ""


		# Create one new volume per mask value per hemisphere, and run atlasquery against each of those volumes:
		for maskValue in ${maskValues}; do
			echo ""
			echo "############################################################"
			echo "# MODEL: ${parentDir}/${model}"
			echo "#     CLUSTER MASK: ${clusterMask}"
			echo "#          MASK VALUE: ${maskValue}"
			echo "############################################################"
			echo ""

			echo "Because FSL's atlasquery doesn't allow user to point"
			echo "at a specific mask intensity of a multi-intensity"
			echo "mask, we must first isolate each mask value into its own"
			echo "file. While doing this we also split the mask into"
			echo "left and right hemispheres because many of"
			echo "atlasquery's a priori regions don't have separate left"
			echo "and right hemisphere compartments:"

			rm -f ${clusterMask}-1mm-?H-maskOrigValue${maskValue}*
			fslmaths ${clusterMask}-1mm -thr ${maskValue} -uthr ${maskValue} -mul ${blockLH} ${clusterMask}-1mm-LH-maskOrigValue${maskValue} -odt char
			fslmaths ${clusterMask}-1mm -thr ${maskValue} -uthr ${maskValue} -mul ${blockRH} ${clusterMask}-1mm-RH-maskOrigValue${maskValue} -odt char
			ls -al ${clusterMask}-1mm-?H-maskOrigValue${maskValue}*

			echo ""
			echo "CREATING LISTS of atlas regions that overlap with ${clusterMask}-1mm-?H-maskOrigValue${maskValue}:"
			for hem in LH RH; do
				echo ""
				echo "############################################################"
				echo "# MODEL: ${parentDir}/${model}"
				echo "#     CLUSTER MASK: ${clusterMask}"
				echo "#          MASK VALUE: ${maskValue}"
				echo "#                HEMISPHERE: ${hem}"
				echo "############################################################"
				echo ""
				echo "~~~~~~~~~~~~~~~ START RESULTS: mask value ${maskValue}, hemisphere ${hem} ~~~~~~~~~~~~~~~"
				atlas=hoCortical
				atlasquery -a "Harvard-Oxford Cortical Structural Atlas"    -m ${clusterMask}-1mm-${hem}-maskOrigValue${maskValue} > ${parentDir}/${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				# create a version of the text file that is prepended with colon-separated identifiers for easy importing into spreadsheets/SPSS/R:
				# (colon-separated because that's the convention that FSL atlasquery uses already)
				rm -f ${parentDir}/colonSeparatedValues-${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				rowIdentifiersForPrepend="${model}:${clusterMask}:${maskValue}:${hem}:${atlas}"
				while read -r line; do
					echo "${rowIdentifiersForPrepend}:${line}"
				done <${parentDir}/${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt > ${parentDir}/colonSeparatedValues-${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt

				atlas=hoSubcortical
				atlasquery -a "Harvard-Oxford Subcortical Structural Atlas" -m ${clusterMask}-1mm-${hem}-maskOrigValue${maskValue} > ${parentDir}/${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				# create a version of the text file that is prepended with colon-separated identifiers for easy importing into spreadsheets/SPSS/R:
				# (colon-separated because that's the convention that FSL atlasquery uses already)
				rm -f ${parentDir}/colonSeparatedValues-${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				rowIdentifiersForPrepend="${model}:${clusterMask}:${maskValue}:${hem}:${atlas}"
				while read -r line; do
					echo "${rowIdentifiersForPrepend}:${line}"
				done <${parentDir}/${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt > ${parentDir}/colonSeparatedValues-${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt

				echo ""	
				echo "DEBUG: Confirm that we have the same number of lines in the prefixed output as the original output:"	
				ls -al  ${parentDir}/hoCortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				wc -l  ${parentDir}/hoCortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				ls -al  ${parentDir}/colonSeparatedValues-hoCortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				wc -l  ${parentDir}/colonSeparatedValues-hoCortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				ls -al  ${parentDir}/hoSubcortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				wc -l  ${parentDir}/hoSubcortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				ls -al  ${parentDir}/colonSeparatedValues-hoSubcortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				wc -l  ${parentDir}/colonSeparatedValues-hoSubcortical-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt
				echo ""

				# generate copy-and-paste commands for manually auditing the most granular version of the output:
				for atlas in hoCortical hoSubcortical; do
					echo ""
					echo "############################################################"
					echo "# MODEL: ${parentDir}/${model}"
					echo "#     CLUSTER MASK: ${clusterMask}"
					echo "#          MASK VALUE: ${maskValue}"
					echo "#                HEMISPHERE: ${hem}"
					echo "#                       ATLAS: ${atlas}"
					echo "############################################################"
					echo ""
					echo "To manually audit overlap list for HEMISPHERE ${hem} in ATLAS ${atlas}"
					echo "just copy-and-paste these commands:"
					echo "      1) cat ${parentDir}/${atlas}-${model}+++${clusterMask}+++${hem}+++cluster${maskValue}.txt"
					echo "      2) fslview -m ortho,lightbox \\"
					echo "         ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -t .4 \\"
					echo "         ${parentDir}/${model}/${clusterMask}-1mm-${hem}-maskOrigValue${maskValue} \\"
					echo "         ${FSLDIR}/data/atlases/HarvardOxford/HarvardOxford-cort-maxprob-thr25-1mm.nii.gz -l MGH-Cortical -t .1 &"
					echo "             (then Tools -> Toolbars -> Atlas Tools)"
					echo "      3) fslview -m ortho,lightbox \\"
					echo "         ${FSLDIR}/data/standard/MNI152_T1_2mm.nii.gz -t .4 \\"
					echo "         ${parentDir}/${model}/${clusterMask} -l \"Random-Rainbow\" &"
					echo ""
				done
				echo "~~~~~~~~~~~~~~~ END RESULTS: mask value ${maskValue}, hemisphere ${hem} ~~~~~~~~~~~~~~~"
			done
		done
	done
done

test

test