By default, mri_glmfit-sim will compute a spatial average inside of a
cluster to create the OCN table, so the values do not agree with the
table below. In he case of area or volume inputs, it needs to compute
the sum instead of the average. I have just modified this program to do
this. I've attached this new version. Copy it to $FREESURFER/bin
To run, add --spatial-sum to the command line
On 1/26/2021 12:56 PM, Na, Xiaoxu wrote:
External Email - Use Caution
To make it clear, I would like to extract the area of each subject in
the cluster. May I know if this is the correct way? If not, it’s
appreciated if you would point it out. Thanks!
*From:* freesurfer-boun...@nmr.mgh.harvard.edu
<freesurfer-boun...@nmr.mgh.harvard.edu> *On Behalf Of *Na, Xiaoxu
*Sent:* Tuesday, January 26, 2021 11:44 AM
*To:* freesurfer@nmr.mgh.harvard.edu
*Subject:* [Freesurfer] Unit of the Numbers in cache.th40.neg.y.ocn.dat
* External Email - Use Caution *
Dear Freesurfer Experts,
As I’d like to extract the average value of each subject in the
cluster after running group analysis, may I know the unit of these
numbers?
Since the one I’m running is surface area, and the cluster is 444mm^2,
all the numbers in the file cache.th40.neg.y.ocn.dat is in the range
of [0.55 1.05], may I know if this is correct and if so, how the unit
is for those? Thanks a lot!
ClusterNo
Max
VtxMax
Size(mm^2)
MNIX
MNIY
MNIZ
CWP
CWPLow
CWPHi
NVtxs
WghtVtx
1
-4.8321
40067
444.11
-58.1
-16.6
-23.9
0.0002
0
0.0004
670
-2784.11
Regards,
Xiaoxu
------------------------------------------------------------------------
Confidentiality Notice: This e-mail message, including any
attachments, is for the sole use of the intended recipient(s) and may
contain confidential and privileged information. Any unauthorized
review, use, disclosure or distribution is prohibited. If you are not
the intended recipient, please contact the sender by reply e-mail and
destroy all copies of the original message.
_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
#!/bin/tcsh -f
set VERSION = 'mri_glmfit-sim @FS_VERSION@';
setenv LANG en_US.UTF-8
set inputargs = ($argv);
set glmdir = ()
set nulltype = ();
set nsim = ()
set thresh = ()
set csdbase = ()
set simsign = ()
set cwpvalthresh = .05
set Overwrite = 0;
set DoSim = 1;
set LF = ();
set AllDoneFile = ();
set Seed = ();
set fwhmOverride = ();
set fwhmAdd = 0;
set UseUniformPDF = 0;
set UniformPDFMin = ();
set UniformPDFMax = ();
set DiagCluster = 0;
set PBNodeType = ();
set UseCache = 0;
set PermForce = 1;
set PermSignFlip = 0; # for osgm-type contrasts
set PermResid = 1; # use residual when running permutation, ter Braak method
set UseGRF = 0;
set volsubject = fsaverage;
set subjectOverride = ();
set Bonferroni = 0;
set ReportCentroid = 0;
set annot = aparc
set IsPVR = 0;
set DoClusterMean = 1
set OutputAnnot = 1
set DoGRFOCNAnat = 0;
set y = ()
set framemask = ()
set PermNonStatCor = 0;
set fwhmmap = ()
set DoSpatialSum = 0
set nJobs = 1;
set DoBackground = 0;
set DoPBSubmit = 0;
set DoPoll = 0;
set tSleepSec = 10;
set CacheDir = $FREESURFER_HOME/average/mult-comp-cor
set CacheLabel = cortex;
set tmpdir = ();
set CleanUp = 1;
set PrintHelp = 0;
if($#argv == 0) goto usage_exit;
set n = `echo $argv | grep -e --help | wc -l`
if($n != 0) then
set PrintHelp = 1;
goto usage_exit;
endif
set n = `echo $argv | grep -e --version | wc -l`
if($n != 0) then
echo $VERSION
exit 0;
endif
source $FREESURFER_HOME/sources.csh
goto parse_args;
parse_args_return:
goto check_params;
check_params_return:
if(! -e $glmdir) then
echo "ERROR: cannot find $glmdir"
exit 1;
endif
set glmfitlog = $glmdir/mri_glmfit.log
if(! -e $glmfitlog) then
echo "ERROR: cannot find $glmfitlog"
exit 1;
endif
if($nulltype != perm) then
set fwhmfile = $glmdir/fwhm.dat
if(! -e $fwhmfile) then
echo "ERROR: cannot find $fwhmfile"
exit 1;
endif
# This may be overridden
set fwhm = `cat $fwhmfile`;
else
set fwhm = 0;
endif
set glmfitcwd = `cat $glmfitlog | awk '{if($1 == "cwd") print $2}'`
if(0 && ! -e $glmfitcwd) then
echo "ERROR: cannot find $glmfitcwd"
exit 1;
endif
set glmfitcwd = `pwd`;
set anattype = volume;
set subject = ();
set hemi = ();
set surf = "white";
set wls = ();
set glmfitcwd = `cat $glmfitlog | awk '{if($1 == "cwd") print $2}'`
# Go through the original mri_glmfit command-line
set glmfitcmd0 = `cat $glmfitlog | awk '{if($1 == "cmdline") print $0}'`
set glmfitcmd = ($glmfitcmd0);
echo $glmfitcmd
set gd2mtx = dods
set UseTable = 0;
set label = ();
while($#glmfitcmd)
set flag = $glmfitcmd[1]; shift glmfitcmd;
#echo $flag
switch($flag)
case "--surf"
case "--surface"
set subject = $glmfitcmd[1]; shift glmfitcmd;
set hemi = $glmfitcmd[1]; shift glmfitcmd;
set anattype = surface;
if($#glmfitcmd != 0) then
set c = `echo $glmfitcmd[1] | cut -c 1-2`;
if("$c" != "--") then
set surf = $glmfitcmd[1]; shift glmfitcmd;
endif
endif
breaksw
case "--fsgd"
shift glmfitcmd;
if($#glmfitcmd > 0) then
if($glmfitcmd[1] == doss) set gd2mtx = doss
endif
breaksw
case "--label"
set label = $glmfitcmd[1]; shift glmfitcmd;
breaksw
case "--perm-force"
# Probably never gets here
set PermForce = 1;
breaksw
case "--table"
case "--y"
if($#y == 0) set y = $glmfitcmd[1];
shift glmfitcmd;
if($DoClusterMean) then
if(! -e $y) then
set y = $glmfitcwd/$y
if(! -e $y) then
set y = $glmdir/`basename $y`
if(! -e $y) then
echo "ERROR: cannot find $y"
echo "If you do not have this file, run with --no-y as the first
option"
exit 1;
endif
endif
endif
endif
if("$flag" == "--table") set UseTable = 1;
breaksw
case "--frame-mask"
set framemask = $glmfitcmd[1];
shift glmfitcmd;
if(! -e $framemask) then
echo "ERROR: cannot find frame mask $framemask"
exit 1;
endif
echo "Using frame mask $framemask"
breaksw
case "--wls"
if($nulltype == perm) then
echo "ERROR: cannot use permutation with weighted least squares (--wls)"
echo "You must re-run mri_glmfit without the --wls flag"
exit 1;
endif
set wls = $glmfitcwd/$glmfitcmd[1]; shift glmfitcmd;
if(! -e $wls) then
set wls = $glmdir/`basename $wls`
if(! -e $wls) then
echo "ERROR: cannot find $wls"
exit 1;
endif
endif
breaksw
case "--osgm"
# Technically don't need to do this here because automatically done in
mri_glmfit
set PermSignFlip = 1;
breaksw
# Ignore these options that have no args
case "cmdline"
case "mri_glmfit"
case "--prune"
case "--no-prune"
case "--pca"
case "--synth"
case "--allowsubjrep"
case "--illcond"
case "--debug"
case "--synth"
case "--cortex"
case "--kurtosis"
case "--nii"
case "--nii.gz"
case "--rescale-x"
case "--no-rescale-x"
case "--fisher"
case "--no-pcc"
case "dods"
case "mri_glmfit.bin"
case "--save-eres"
case "--eres-save"
case "--save-fwhm-map"
case "--mgz"
case "--mgh"
case "--nii"
case "--nii.gz"
breaksw
# Ignore these options that have one arg
case "--glmdir"
case "--o"
case "--C"
case "--X"
case "--w"
case "--fwhm"
case "--mask"
case "--seed"
case "--var-fwhm"
case "--yffxvar"
case "--ffxdof"
case "--ffxdofdat"
case "--exclude-frame"
shift glmfitcmd;
breaksw
case "--pvr"
set IsPVR = 1;
shift glmfitcmd;
breaksw
# Ignore these options that have two args
case "--rand-split"
shift glmfitcmd;
shift glmfitcmd;
breaksw
default
echo ""
echo "INFO: unrecognized mri_glmfit cmd option $flag"
echo " ... probably not a problem"
echo ""
#sleep 10
#exit 1;
breaksw
endsw
end
set mask = ();
if(! $UseTable) then
set mask = `stem2fname $glmdir/mask`
if($status) then
echo "$mask"
exit 1;
endif
endif
if(0 && $UseGRF && $anattype != volume) then
# not true anymore
echo "ERROR: can only use --grf with volume"
exit 1;
endif
if($UseGRF) then
set csdbase = "grf.th$thresh.$simsign"
#set tmp = `echo "$thresh*10"|bc -l`
#set tmp = `printf %02d $tmp`
#set csdbase = "grf.th$tmp.$simsign"
endif
if($#subjectOverride) set subject = $subjectOverride
# Make sure the subject exists
if($#subject != 0) then
if(! -e $SUBJECTS_DIR/$subject) then
echo "ERROR: cannot find $subject in $SUBJECTS_DIR"
exit 1;
endif
echo "SURFACE: $subject $hemi"
endif
# Get a list of contrasts, etc
set FSGDhasCon = 0;
if(-e $glmdir/y.fsgd) then
# check whether contrasts were specified in the FSGD
set FSGDhasCon = `cat $glmdir/y.fsgd | awk '{if($1 == "Contrast") print 1}' |
wc -l`
endif
if($#tmpdir == 0) set tmpdir = $glmdir/tmp.mri_glmfit-sim-$$
mkdir -p $tmpdir
set clist = ($glmdir/*/C.dat)
if($status) then
echo "ERROR: cannot find any contrasts"
exit 1;
endif
set clist2 = ();
set conlist = ();
foreach c ($clist)
set tmp = `dirname $c`
set conname = `basename $tmp`
set conlist = ($conlist $conname);
set confile = $tmpdir/$conname.mtx
cp $c $confile
# If Cons in the FSGD, dont spec on glmfit cmd line
if($FSGDhasCon == 0) set clist2 = ($clist2 --C $confile);
# Make sure sig is there
set sig = `stem2fname $glmdir/$conname/sig`
if($status) then
echo "$sig"
exit 1;
endif
# Make sure that csd is not there
if($DoSim) then
ls $glmdir/csd/$csdbase.j???-$conname.csd >& /dev/null
set csdexists = (! $status)
if($csdexists) then
set csdfiles = ($glmdir/csd/$csdbase.j???-$conname.csd)
if(! $Overwrite) then
echo "ERROR: $csdfiles already exists"
echo "Delete it or run with --overwrite"
exit 1;
endif
rm -f $csdfiles
endif
endif
end
set StartDate = `date`
# Create log file
if($#LF == 0) then
if($UseCache == 0) then
set LF = $glmdir/$csdbase.mri_glmfit-sim.log
else
set LF = $glmdir/cache.mri_glmfit-sim.log
endif
rm -f $LF
endif
echo "log file is $LF"
echo "" | tee -a $LF
echo "cd `pwd`" | tee -a $LF
echo $0 | tee -a $LF
echo $inputargs | tee -a $LF
echo "" | tee -a $LF
echo $VERSION | tee -a $LF
date | tee -a $LF
uname -a | tee -a $LF
echo $user | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "FREESURFER_HOME $FREESURFER_HOME" | tee -a $LF
echo "" | tee -a $LF
echo "Original mri_glmfit command line:" | tee -a $LF
echo $glmfitcmd0 | tee -a $LF
echo "" | tee -a $LF
echo "DoSim = $DoSim" | tee -a $LF
echo "UseCache = $UseCache" | tee -a $LF
echo "DoPoll = $DoPoll" | tee -a $LF
echo "DoPBSubmit = $DoPBSubmit" | tee -a $LF
echo "DoBackground = $DoBackground" | tee -a $LF
echo "DiagCluster = $DiagCluster" | tee -a $LF
echo "gd2mtx = $gd2mtx" | tee -a $LF
if($#Seed) echo "Seed = $Seed" | tee -a $LF
echo "fwhm = $fwhm" | tee -a $LF
if($#fwhmOverride) then
# Change fwhm from value in glmdir to override value
echo "fwhmOverride = $fwhmOverride" | tee -a $LF
set fwhm = $fwhmOverride
endif
mkdir -p $glmdir/csd
if($UseCache) then
#set fwhm = `cat $glmdir/fwhm.dat`;
set fwhmStr = `perl -e "printf('"'%02d'"',int ( $fwhm+0.5+$fwhmAdd ) ) "`
if($fwhmStr == 00) set fwhmStr = 01;
set threshStr = `perl -e "print 10*$thresh"`;
set csdCache =
$CacheDir/$subject/$hemi/$CacheLabel/fwhm$fwhmStr/$simsign/th$threshStr/mc-z.csd
if(! -e $csdCache) then
echo "ERROR: cannot find $csdCache" | tee -a $LF
echo " This often happens when the measured FWHM is very high and their is
not a table entry." | tee -a $LF
echo " Consider using lower smoothing levels." | tee -a $LF
exit 1;
endif
echo CSD $csdCache | tee -a $LF
if($#csdbase == 0) set csdbase = cache.th$threshStr.$simsign
endif
if($DoSim) then
# Use fsgd if there, if not use X
set fsgd = ();
set X = ();
if(-e $glmdir/y.fsgd) then
set fsgd = $glmdir/y.fsgd
else
set X = $glmdir/Xg.dat
endif
# Number of sim iterations per job
set nSimPerJob = `echo "$nsim/$nJobs" | bc`;
echo "nSimPerJob = $nSimPerJob"| tee -a $LF
if($DoPoll) then
set PollDir = $glmdir/csd/poll
mkdir -p $PollDir
# User can create this file to stop simulation
set StopFile = $PollDir/StopSim
rm -f $StopFile
endif
# ----------------- Run glmfit simulation ----------------------
# This may take a while
@ nthJob = 1
while($nthJob <= $nJobs)
set jpad = `printf j%03d $nthJob`
echo "$nthJob/$nJobs `date`" | tee -a $LF
set cmd = (mri_glmfit $clist2)
set cmd = ($cmd --sim $nulltype $nSimPerJob $thresh \
$glmdir/csd/$csdbase.$jpad)
set cmd = ($cmd --sim-sign $simsign --fwhm $fwhm)
if($#fsgd) set cmd = ($cmd --fsgd $fsgd $gd2mtx)
if($#X) set cmd = ($cmd --X $X)
if($#wls) set cmd = ($cmd --wls $wls)
if($#framemask) set cmd = ($cmd --frame-mask $framemask)
# --y must go after --fsgd
if(! $UseTable) then
set cmd = ($cmd --mask $mask)
if($PermResid == 0) set cmd = ($cmd --y $y)
if($PermResid == 1) set cmd = ($cmd --y $residual)
if($PermSignFlip == 1) set cmd = ($cmd --perm-1)
if($PermNonStatCor) set cmd = ($cmd --perm-nonstatcor)
endif
if($UseTable) set cmd = ($cmd --table $y)
#if($#label) set cmd = ($cmd --label $label)
if($PermForce) set cmd = ($cmd --perm-force)
if($anattype == surface) set cmd = ($cmd --surf $subject $hemi $surf);
if($#Seed) set cmd = ($cmd --seed $Seed);
if($UseUniformPDF) set cmd = ($cmd --uniform $UniformPDFMin $UniformPDFMax);
if($DiagCluster) set cmd = ($cmd --diag-cluster --debug)
if($?verbose) set cmd = ($cmd --debug)
if($DoPoll) then
set DoneFile = $PollDir/done.$csdbase.$jpad
rm -f $DoneFile
set cmd = ($cmd --sim-done $DoneFile)
endif
echo $cmd | tee -a $LF
if(! $DoPoll) then
# Run directly in shell
$cmd | tee -a $LF
if($status && ! $DiagCluster) exit 1;
date | tee -a $LF
# ----------------------------------------------------
if($DiagCluster) then
set conname = $conlist[1];
set sig = ./$conname-sig.mgh
set csd = $glmdir/csd/$csdbase.$jpad-$conname.csd
set sum = $glmdir/csd/$csdbase.$jpad-$conname.diag
if($anattype == surface) then
set cmd = (mri_surfcluster --subject $subject \
--hemi $hemi --surf $surf --cwpvalthresh $cwpvalthresh)
if($ReportCentroid) set cmd = ($cmd --centroid);
set c = 4;
endif
if($anattype == volume) then
set cmd = (mri_volcluster)
set c = 3;
endif
set cmd = ($cmd --in $sig --mask $mask --cwpvalthresh $cwpvalthresh \
--sum $sum --thmin $thresh --sign $simsign)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
set maxclustB = `grep -v \# $csd | awk '{print $3}' | sort -nr | head
-n 1`;
set maxclustA = `grep -v \# $sum | awk -v c=$c '{print $c}' | sort -nr
| head -n 1`;
echo "MaxClusterSizeA $maxclustA"| tee -a $LF
echo "MaxClusterSizeB $maxclustB"| tee -a $LF
exit 1;
endif
# ----------------------------------------------------
else
if($DoBackground) then
# Run in background
$cmd >>& $LF &
endif
if($DoPBSubmit) then
# pbsubmit
echo "PBSUBMIT ---------------------------------" | tee -a $LF
pbsubmit $PBNodeType -c "$cmd >> $LF" | tee -a $LF
sleep 10;
echo "" | tee -a $LF
endif
endif
@ nthJob = $nthJob + 1
end
if($DoPoll) then
# Possible bug: if job above dies without creating done file
echo "Polling" | tee -a $LF
@ nthPoll = 0;
set Done = 0;
while(! $Done)
if(-e $StopFile) then
# Something has gone wrong, maybe user deleted to halt
echo "WARNING: simulation aborted by user" | tee -a $LF
exit 1;
endif
if(! -e $PollDir) then
# Something has gone wrong, maybe user deleted to halt
echo "ERROR: cannot find $PollDir, aborting" | tee -a $LF
exit 1;
endif
@ nthPoll = $nthPoll + 1;
set Done = 1;
@ nthJob = 1
while($nthJob <= $nJobs)
set jpad = `printf j%03d $nthJob`
set DoneFile = $PollDir/done.$csdbase.$jpad
if(! -e $DoneFile) then
echo "Poll $nthPoll job $nthJob `date`" | tee -a $LF
set Done = 0;
break;
endif
@ nthJob = $nthJob + 1
end
sleep $tSleepSec
end
endif
endif #DoSim
if($DoPBSubmit) then
# pbsubmit the surfcluster part
set cmd = (mri_glmfit-sim --glmdir $glmdir --tmp $tmpdir)
set cmd = ($cmd --sim-sign $simsign --no-sim $csdbase --cwpvalthresh
$cwpvalthresh)
set cmd = ($cmd --log $LF)
if($PermNonStatCor) set cmd = ($cmd --perm-nonstatcor)
echo "------------------------------------------------" | tee -a $LF
echo $cmd | tee -a $LF
pbsubmit $PBNodeType -c "$cmd"
sleep 1;
echo ""
echo "Just submitted this as a job, and exiting now."
echo ""
# Could wait until it is done by polling for AllDoneFile
exit 0;
endif
# -------------- now run clustering -------------------- #
foreach conname ($conlist)
if($UseTable) continue;
set csdlist = ();
if(! $UseGRF) then
if($UseCache) then
set csdfiles = ($csdCache)
else
ls $glmdir/csd/$csdbase.j???-$conname.csd >& /dev/null
if($status) then
echo "ERROR: cannot find any csd files"| tee -a $LF
exit 1;
endif
set csdfiles = ($glmdir/csd/$csdbase.j???-$conname.csd)
endif
set csdlist = ();
foreach csd ($csdfiles)
set csdlist = ($csdlist --csd $csd);
end
endif
set sig = `stem2fname $glmdir/$conname/sig`
set ext = `fname2ext $sig`
set vwsig = $glmdir/$conname/$csdbase.sig.voxel.$ext
set vwsigmax = $glmdir/$conname/$csdbase.sig.voxel.max.dat
set cwsig = $glmdir/$conname/$csdbase.sig.cluster.$ext
set msig = $glmdir/$conname/$csdbase.sig.masked.$ext
set ocn = $glmdir/$conname/$csdbase.sig.ocn.$ext
set oannot = $glmdir/$conname/$csdbase.sig.ocn.annot
set sum = $glmdir/$conname/$csdbase.sig.cluster.summary
set csdpdf = $glmdir/$conname/$csdbase.pdf.dat
if($anattype == surface) then
set cmd = (mri_surfcluster --in $sig --mask $mask \
--cwsig $cwsig --sum $sum --ocn $ocn \
--annot $annot --cwpvalthresh $cwpvalthresh --o $msig --no-fixmni)
set cmd = ($cmd --csd-out $glmdir/csd/all.$csdbase-$conname.csd)
if(! $UseGRF) set cmd = ($cmd $csdlist --csdpdf $csdpdf --vwsig $vwsig
--vwsigmax $vwsigmax)
if($UseGRF) set cmd = ($cmd --fwhm $fwhm --hemi $hemi --subject $subject
--thmin $thresh --sign $simsign)
if($PermNonStatCor) set cmd = ($cmd --fwhm-map $fwhmmap)
if($OutputAnnot) set cmd = ($cmd --oannot $oannot)
if($Bonferroni) set cmd = ($cmd --bonferroni $Bonferroni)
if($ReportCentroid) set cmd = ($cmd --centroid);
set cmd = ($cmd --surf $surf); # bug to have to do this
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
if($DoClusterMean) then
# Create an average "waveform" for each cluster
set avgwf = $glmdir/$conname/$csdbase.y.ocn.dat
set cmd = (mri_segstats --seg $ocn --exclude 0 \
--i $y --sum /tmp/mri_glmfit-sim.junk.$$)
if($DoSpatialSum == 0) set cmd = ($cmd --avgwf $avgwf)
if($DoSpatialSum == 1) set cmd = ($cmd --sumwf $avgwf)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
rm -f /tmp/mri_glmfit-sim.junk.$$
endif
endif
if($anattype == volume) then
mri_info --o /tmp/mri_glmfit-sim.junk.$$.cres --cres $sig
set voxsize = `cat /tmp/mri_glmfit-sim.junk.$$.cres`
if($voxsize != 2) then
echo "ERROR: voxel size of input is $voxsize, expecting 2"
exit 1
endif
set reg = $FREESURFER_HOME/subjects/$volsubject/mri.2mm/reg.2mm.dat
set aseg = aparc+aseg.mgz
if(! -e $FREESURFER_HOME/subjects/$volsubject/mri/$aseg) set aseg = aseg.mgz
set cmd = (mri_volcluster --in $sig --mask $mask --reg $reg --no-fixmni\
--cwsig $cwsig --sum $sum --ocn $ocn --cwpvalthresh $cwpvalthresh \
--seg $volsubject $aseg --out $msig --allowdiag)
if($Bonferroni) set cmd = ($cmd --bonferroni $Bonferroni)
if(! $UseGRF) set cmd = ($cmd --csdpdf $csdpdf $csdlist --vwsig $vwsig
--vwsigmax $vwsigmax)
if($UseGRF) set cmd = ($cmd --fwhmdat $glmdir/fwhm.dat --sign $simsign
--thmin $thresh --sign $simsign)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
if($DoClusterMean) then
# Create an average "waveform" for each cluster
set avgwf = $glmdir/$conname/$csdbase.y.ocn.dat
set cmd = (mri_segstats --seg $ocn --exclude 0 \
--i $y --sum /tmp/mri_glmfit-sim.junk.$$)
if($DoSpatialSum == 0) set cmd = ($cmd --avgwf $avgwf)
if($DoSpatialSum == 1) set cmd = ($cmd --sumwf $avgwf)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
rm -f /tmp/mri_glmfit-sim.junk.$$
endif
if($UseGRF && $DoGRFOCNAnat) then
# Map the OCN into anat space, takes a lot of disk space
set ocnanat = $glmdir/$conname/$csdbase.sig.ocn.anat.$ext
set cmd = (mri_vol2vol --interp nearest --mov $ocn --regheader \
--targ $SUBJECTS_DIR/fsaverage/mri/orig.mgz \
--o $ocnanat --no-save-reg)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
endif
endif
# For some reason cwsig has 4 frames, just take 1st
set cmd = (mri_convert $cwsig $cwsig --frame 0)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
end
if($CleanUp) rm -r $tmpdir
if($#AllDoneFile) touch $AllDoneFile
echo $StartDate | tee -a $LF
date | tee -a $LF
echo "mri_glmfit-sim done" | tee -a $LF
exit 0
###############################################
############--------------##################
parse_args:
set cmdline = ($argv);
while( $#argv != 0 )
set flag = $argv[1]; shift;
switch($flag)
case "--glmdir":
if( $#argv < 0) goto arg1err;
set glmdir = $argv[1]; shift;
breaksw
case "--y":
if($#argv < 0) goto arg1err;
# Override input in fsgd file
set y = $argv[1]; shift;
if(! -e $y) then
echo "ERROR: cannot find $y"
exit 1
endif
breaksw
case "--mczsim":
case "--mcsim":
case "--cache":
if( $#argv < 2) goto arg2err;
set thresh = $argv[1]; shift;
set simsign = $argv[1]; shift;
set UseCache = 1;
set DoSim = 0;
# Handles the case where LC_NUMERIC is not set to en_US.UTF-8
set thresh = `echo $thresh | sed 's/,/./g'`
set thresh = `printf %2.1f $thresh`
if($thresh != 1.3 && $thresh != 2.0 && $thresh != 2.3 && \
$thresh != 3.0 && $thresh != 3.3 && $thresh != 4.0) then
echo "ERROR: thresh = $thresh, must be 1.3, 2.0, 2.3, 3.0, 3.3, 4.0"
exit 1;
endif
breaksw
case "--mczsim-dir":
case "--cache-dir":
if( $#argv < 1) goto arg1err;
set CacheDir = $argv[1]; shift;
set UseCache = 1;
set DoSim = 0;
breaksw
case "--mczsim-label":
case "--cache-label":
if( $#argv < 1) goto arg1err;
set CacheLabel = $argv[1]; shift;
set UseCache = 1;
set DoSim = 0;
breaksw
case "--vol-subject":
if( $#argv < 1) goto arg1err;
set volsubject = $argv[1]; shift;
breaksw
case "--grf":
if( $#argv < 2) goto arg2err;
set thresh = $argv[1]; shift;
set simsign = $argv[1]; shift;
#if($simsign != pos && $simsign != neg) then
#echo "ERROR: must use pos or neg with --grf"
#exit 1;
#endif
set UseGRF = 1;
set DoSim = 0;
breaksw
case "--perm":
if( $#argv < 3) goto arg4err;
set nulltype = perm;
set nsim = $argv[1]; shift;
set thresh = $argv[1]; shift;
set simsign = $argv[1]; shift;
set threshStr = `perl -e "print 10*$thresh"`;
set csdbase = perm.th$threshStr.$simsign
breaksw
case "--perm-nonstatcor":
set PermNonStatCor = 1;
set csdbase = perm.nsc.th$threshStr.$simsign
breaksw
case "--sim":
if( $#argv < 4) goto arg4err;
set nulltype = $argv[1]; shift;
set nsim = $argv[1]; shift;
set thresh = $argv[1]; shift;
set csdbase = $argv[1]; shift;
breaksw
case "--sim-sign":
if( $#argv < 1) goto arg1err;
set simsign = $argv[1]; shift;
breaksw
case "--subject-override":
if( $#argv < 1) goto arg1err;
set subjectOverride = $argv[1]; shift;
breaksw
case "--centroid":
set ReportCentroid = 1;
breaksw
case "--log":
if( $#argv < 1) goto arg1err;
set LF = $argv[1]; shift;
breaksw
case "--a2009s":
set annot = "aparc.a2009s"
breaksw
case "--annot":
if( $#argv < 1) goto arg1err;
set annot = $argv[1]; shift;
breaksw
case "--overwrite":
set Overwrite = 1;
breaksw
case "--pbsubmit":
if($#argv < 1) goto arg1err;
set nJobs = $argv[1]; shift;
set DoPBSubmit = 1;
set DoPoll = 1;
set DoBackground = 0;
breaksw
case "--bg":
if($#argv < 1) goto arg1err;
set nJobs = $argv[1]; shift;
set DoBackground = 1;
set DoPoll = 1;
set DoPBSubmit = 0;
breaksw
case "--no-bg":
set DoBackground = 0;
set DoPoll = 0;
breaksw
case "--sleep":
if($#argv < 1) goto arg1err;
set tSleepSec = $argv[1]; shift;
breaksw
case "--seed":
if($#argv < 1) goto arg1err;
set Seed = $argv[1]; shift;
breaksw
case "--fwhm-override":
if($#argv < 1) goto arg1err;
set fwhmOverride = $argv[1]; shift;
breaksw
case "--fwhm-add":
if($#argv < 1) goto arg1err;
set fwhmAdd = $argv[1]; shift;
breaksw
case "--perm-force":
set PermForce = 1;
breaksw
case "--perm-resid":
set PermResid = 1;
set PermForce = 1; # purpose of using residual is to allow non-orthog
breaksw
case "--no-perm-force":
set PermForce = 0;
breaksw
case "--no-perm-resid":
set PermForce = 0;
set PermResid = 0;
breaksw
case "--tmp":
if($#argv < 1) goto arg1err;
set tmpdir = $argv[1]; shift;
breaksw
case "--bonferroni":
if($#argv < 1) goto arg1err;
set Bonferroni = $argv[1]; shift;
breaksw
case "--2spaces":
set Bonferroni = 2;
breaksw
case "--3spaces":
set Bonferroni = 3;
breaksw
case "--cwpvalthresh":
case "--cwp":
if($#argv < 1) goto arg1err;
set cwpvalthresh = $argv[1]; shift;
breaksw
case "--uniform":
if($#argv < 2) goto arg2err;
set UniformPDFMin = $argv[1]; shift;
set UniformPDFMax = $argv[1]; shift;
set UseUniformPDF = 1;
breaksw
case "--base":
if($#argv < 1) goto arg1err;
set csdbase = $argv[1]; shift;
breaksw
case "--no-sim":
if($#argv < 1) goto arg1err;
set csdbase = $argv[1]; shift;
set DoSim = 0;
breaksw
case "--spatial-sum":
set DoSpatialSum = 1;
breaksw
case "--diag-cluster":
set DiagCluster = 1;
breaksw
case "--no-out-annot":
set OutputAnnot = 0
breaksw
case "--no-cluster-mean":
case "--no-y":
set DoClusterMean = 0;
breaksw
case "--no-grf-ocn-anat":
set DoGRFOCNAnat = 0;
breaksw
case "--grf-ocn-anat":
set DoGRFOCNAnat = 1;
breaksw
case "--debug":
set verbose = 1;
set echo = 1;
breaksw
default:
echo ERROR: Flag $flag unrecognized.
echo $cmdline
exit 1
breaksw
endsw
end
goto parse_args_return;
############--------------##################
############--------------##################
check_params:
if($#glmdir == 0) then
echo "ERROR: must spec --glmdir"
exit 1;
endif
if($DoSim) then
if($#nulltype == 0) then
echo "ERROR: must spec --sim"
exit 1;
endif
if($#simsign == 0) then
echo "ERROR: must spec --sim-sign"
exit 1;
endif
endif
set residual = ()
if($nulltype == perm && $PermResid) then
set stem = $glmdir/eres
set residual = `stem2fname $stem`
if($status) then
echo "ERROR: cannot find residual $stem needed with --perm-resid"
echo " when running mri_glmfit, make sure to include --eres-save"
exit 1;
endif
endif
if($DiagCluster && $DoBackground) then
echo "ERROR: cannot background with --diag-cluster"
exit 1;
endif
if($DiagCluster && $DoPBSubmit) then
echo "ERROR: cannot pbsubmit with --diag-cluster"
exit 1;
endif
if($DiagCluster && $DoSim == 0) then
echo "ERROR: must simulate with --diag-cluster"
exit 1;
endif
if($PermNonStatCor) then
set stem = $glmdir/fwhm
set fwhmmap = `stem2fname $stem`
if($status) then
echo "ERROR: cannot find fwhmmap $stem needed with --perm-nonstatcor"
echo " when running mri_glmfit, make sure to include --save-fwhm-map"
exit 1;
endif
endif
goto check_params_return;
############--------------##################
############--------------##################
arg1err:
echo "ERROR: flag $flag requires one argument"
exit 1
############--------------##################
############--------------##################
arg2err:
echo "ERROR: flag $flag requires two arguments"
exit 1
############--------------##################
############--------------##################
arg4err:
echo "ERROR: flag $flag requires four arguments"
exit 1
############--------------##################
############--------------##################
usage_exit:
echo ""
echo "mri_glmfit-sim --help"
echo ""
echo " --glmdir glmdir"
echo " --cwp cwthresh : cluster-wise p-value threshold ($cwpvalthresh)"
echo ""
echo "Use pre-computed z-based monte-carlo simulations (see mri_mcsim)"
echo " --mczsim vwthreshold sign"
echo "If you have run your own simulations (mri_mcsim), then specify "
echo " --mczsim-dir dir : default is FREESURFER_HOME/average/mult-comp-cor"
echo " --mczsim-label label : default is cortex"
echo ""
echo "Use permutation simulation"
echo " --perm nsim CFT sign : permutation simulation with nsim iterations"
echo " --perm-resid : add this with non-orthogonal designs to run
permutation"
echo " on the residual rather than on the raw data (must have run
mri_glmfit"
echo " with --eres-save)"
echo " --perm-signflip : use sign flipping instead of shuffling (automatic
for --osgm)"
echo ""
echo "Use Gaussian Random Fields (GRF), volumes only"
echo " --grf vwthreshold sign : -log10(p) and sign (pos, neg, abs)"
echo ""
echo "Correct for multiple spaces"
echo " --2spaces : additional Bonferroni correction across 2 spaces (eg, lh,
rh)"
echo " --3spaces : additional Bonferroni correction across 3 spaces (eg, lh,
rh, mni305)"
echo ""
echo " --overwrite : delete previous CSDs"
echo " --bg njobs : divide sim into njobs and put in background"
echo " --sleep tSleepSec : number of seconds to sleep between background
polls "
echo ""
echo " --a2009s : use aparc.a2009s instead of aparc for reporting region of
vertex max"
echo " --annot annot : use annot instead of aparc for reporting region of
vertex max"
echo ""
echo " --log logfile : default is csdbase.mri_glmfit-sim.log"
echo ""
echo " --base csdbase : override csdbase name"
echo " --no-sim csdbase : do not simulate, only run cluster"
echo " --seed seed : set simluation seed"
echo " --fwhm-override fwhm : override fwhm in glmdir":
echo " --fwhm-add fwhmAdd : add fwhmAdd to the estimated fwhm before
addressing table":
echo " --uniform min max : use uniform PDF instead of gaussian":
echo " --no-out-annot : do not output a cluster annotation"
echo " --no-cluster-mean : do not compute means of each subject in each
cluster (was --no-y)"
echo " --y yFile : use yFile instead of determining yFile from glmdir"
echo " --centroid : report the coordinates/annotation of the centroid instead
of max"
echo " --spatial-sum : compute the sum over voxels in the cluster rather than"
echo " the average when creating the y.ocn file. This is good for when the
input is area or volume"
echo " --help"
echo ""
if(! $PrintHelp) exit 1;
echo $VERSION
cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'
exit 1;
#---- Everything below here is printed out as part of help -----#
BEGINHELP
This program is a script to run corrections for multiple comparisons
(CMC) on the volume or the surface. It is a front-end for the
mri_glmfit, mri_surfcluster, and mri_volcluster programs.
To use this program:
1. Run mri_glmfit on your data, giving it the glm directory
(--glmdir), and all the other arguments.
2. Choose a voxel-wise (cluster-forming) threshold vwthresh. The
vwthresh takes the value -log10(p) where p is the p-value. So if
vwthresh = 2, then p<.01. You should understand the difference
between the cluster-forming threshold and the cluster-wise
threshold (--cwp). For mczsim, vwthresh must be 1.3, 2.0, 2.3, 3.0,
3.3, or 4.0 (unless you have simulated more values). There are no
good methods to choose this value. Typically 2 or 3 is used.
It is not recommended to use 1.3 unless performing a permutation
operation.
3. Choose the sign of the effect, either positive (pos), negative
(neg) or absolute/unsigned (abs). If you have an a priori
hypothesis (eg, group1 > group2), then choose positive. If you do
not expect a particular sign, then use abs. Choosing a sign will
yield stronger results (ie, smaller p-values), but you will lose
voxels with the opposite sign. NOTE: the "abs" for GRF is only
approximate. For volumes, it may give very conservative results.
4. Choose the method. This will be either: (1) --mczsim z-based monte
carlo simulation (See
surfer.nmr.mgh.harvard.edu/fswiki/BuildYourOwnMonteCarlo).
(2) --perm for permutation simulation, or (3) --grf for Gaussian Random Field
method (analytical). For permution, you must choose the number of simulation
iterations (typically 5000-1000); permutation may take a few hours to
complete. GRF is for volumes only.
Now run the command, eg,
mri_glmfit-sim --glmdir glmdir --mczsim 2 abs
or
mri_glmfit-sim --glmdir glmdir --sim perm 1000 3 abs
The output will be in each contrast folder and will be something like
the following (where csdbase will be "cache" with --mczsim or "perm"
with --perm or "grf" with --grf)
csdbase.sig.voxel.mgh
csdbase.sig.cluster.mgh
csdbase.sig.cluster.summary
csdbase.y.ocn.dat
csdbase.sig.ocn.mgh
csdbase.sig.ocn.annot (for surfaces only)
csdbase.sig.masked.mgh
csdbase.sig.voxel.mgh - the sig volume corrected for multiple
comparisons on a voxel-wise basis. The threshhold and sign are
irrelevant. The value at each voxel is the corrected -log10(p-value)
for that voxel.
csdbase.sig.cluster.mgh - the sig volume corrected for multiple
comparisons on a cluster-wise basis. The value at each voxel
is the -log10(p), where p is the pvalue of the cluster at
that voxel. If that voxel does not belong to a cluster, its
value will be 0.
csdbase.sig.cluster.summary - this is the cluster summary table
(a simple text file).
csdbase.y.ocn.dat - this is a summary of the input (y) over each
cluster. It has a column for each cluster. Each row is a subject. The
value is the average of the input (y) in each cluster. This is a
simple text file.
csdbase.sig.ocn.mgh - output cluster number volume. The value of the
voxel is the integer cluster number that that voxel belongs to. This
can be used like a segmentation (eg, to load into tkmedit or
use for running mri_segstats).
csdbase.sig.ocn.annot (for surfaces only) - this is an annotation
with the annotation name being the cluster number. This can be
used like any other annotation (eg, to load into tksurfer or
as input for mri_segstats or mris_anatomical_stats).
csdbase.sig.masked.mgh - the original sig volume masked to show
only clusters.
OTHER OPTIONS
--overwrite
Overwrite existing CSD files when performing simulation. Default
is to not overwrite so that you do not clobber CSD files that you
spent a weekend creating.
--bg njobs
CAREFUL!!!! WARNING!!! Divide the number of simulation interations
into njobs background jobs. If the simulation jobs die uncleanly for
some reason, this will never quit because it looks for a file created
by the simulation job before continuing. If this file is not created,
then it will keep polling for this file forever. ALSO: if you are not
running on a computer with multiple nodes, then this will not make it
run any faster (in fact, it will run slower). You can stop
backgrounded jobs creating a file called glmdir/csd/poll/StopSim. If
you are in the Martinos center, do not attempt to use --bg on the
launchpad computational cluster.
--cwp cwpvalthresh
--cwpvalthresh cwpvalthresh
Only report clusters that have p < cwpvalthresh. Default is .05. If
you want to change this after running a permutation simulation, you do
not need to re-run the entire simulation. Just run mri_glmfit with the
--no-sim option.
--2spaces
--3spaces
This adds an addition correction for multiple comparisons for when the
current correction is part of a larger analysis. For example, if you
are doing a structural study looking at both the left and the right
hemispheres, mri_glmfit-sim must be run separately for each
hemisphere. By default, it will give you p-values corrected only for
one hemisphere. To correct for both hemispheres use --2spaces. For an
fMRI study in which you are using both hemispheres and subcortical
structures, then use --3spaces. The correction applied is Bonferroni
(N=2 for --2spaces and N=3 for --3spaces).
--no-sim csdbase
Run this if you have already run the permutation stimulation but just
want to run the correction for multiple comparisons again to change
the clusterwise threshold. "csdbase" should be "perm"
--seed seed
Set the seed for simulation random number generator. This is mostly
used for debugging.
--fwhm-override fwhm
Override fwhm in glmdir with given value. This is mostly
used for debugging.
--no-cluster-mean
Do not compute means of each subject in each cluster. This computation
requires the y input file which might not be available in all cases.
When using this option, put it first on the command-line.
--y yFile
Specify the GLM input (ie, argument to --y in mri_glmfit) rather
than determining it from the glmdir. This can be handy if the
y file was deleted at some point.
--perm nsim CFT sign
Perform a permutation simulation to correct for multiple comparisons using
nsim iterations, a cluster forming threshold of CFT (in -log10(p) units),
using the given contrast sign (pos, neg, abs).
--perm-signflip
Perform sign flipping instead of shuffling. This is required (and
automatic) when --osgm is used. But it can be helpful when looking
at any main effect. Eg, if you have two groups, then it is usually
the case that you will look at the difference between the groups.
If you want to look at the mean of the groups, then you would
want to enable sign flipping.
--perm-resid
--no-perm-resid
perm-resid is used for non-orthogonal designs (eg, when there is a
covariate present). It runs the permutation using the residual of the
analysis as the input rather than the raw data (the residual will have
the covariate regressed out). You must have originally run mri_glmfit
with the --eres-save option. This method is an APPROXIMATION known as
the ter Braak method (ter Braak, 1992) that allows for permutation to
be run on non-orthogonal designs. See also Winkler, 2014. This is
turned ON by default. Use --no-perm-resid to turn off.
ter Braak, C.J.F., 1992. Permutation versus bootstrap significance
tests in multiple regression and ANOVA. In: Jöckel, K.-H., Rothe, G.,
Sendler, W. (Eds.), Bootstrapping and Related
Techniques. No. 1989. Springer-Verlag, Berlin, pp. 79–86.
Winkler, et al, 2014. Permutation inference for the general linear
model. Neuroimage 92 (2014) 381–397.
_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
