That was an oversight. I've attached a version that should work. Can you
test it for me? Once I know it works I'll check it in to tree.
doug
On 09/03/2014 12:30 PM, Eryilmaz, Huseyin Hamdi wrote:
Dear FS experts,
I have a question about a MC output file. When I run the MC sim
(using mri_glmfit-sim) on the surface, it generates a bunch of files
including a 'sig.masked.mgh', which is nice for displaying the
significance map limited to the clusters that survive the MC
correction. However, when I run it on the volume, it doesn't generate
this specific file. Does anyone know the reason for that? Is there an
alternative way to display the significance map that show only the
clusters that survive the MC correction?
Many thanks!
Hamdi
--
Hamdi Eryilmaz, PhD
Massachusetts General Hospital
A.A. Martinos Center for Biomedical Imaging
Psychiatric Neuroimaging Division
149 13th St, Charlestown, MA 02129
Phone: +1 617 643 7462
Email:heryil...@partners.org <mailto:heryil...@partners.org>
_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
--
Douglas N. Greve, Ph.D.
MGH-NMR Center
gr...@nmr.mgh.harvard.edu
Phone Number: 617-724-2358
Fax: 617-726-7422
Bugs: surfer.nmr.mgh.harvard.edu/fswiki/BugReporting
FileDrop: https://gate.nmr.mgh.harvard.edu/filedrop2
www.nmr.mgh.harvard.edu/facility/filedrop/index.html
Outgoing: ftp://surfer.nmr.mgh.harvard.edu/transfer/outgoing/flat/greve/
#!/bin/tcsh -f
set VERSION = '$Id: mri_glmfit-sim,v 1.36.2.6 2013/06/05 17:06:31 greve Exp $';
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 UseUniformPDF = 0;
set UniformPDFMin = ();
set UniformPDFMax = ();
set DiagCluster = 0;
set PBNodeType = ();
set UseCache = 0;
set PermForce = 0;
set UseGRF = 0;
set volsubject = fsaverage;
set subjectOverride = ();
set Bonferroni = 0;
set ReportCentroid = 0;
set annot = aparc
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
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
set mask = `stem2fname $glmdir/mask`
if($status) then
echo "$mask"
exit 1;
endif
if($nulltype != perm) then
set fwhmfile = $glmdir/fwhm.dat
if(! -e $fwhmfile) then
echo "ERROR: cannot find $fwhm"
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(! -e $glmfitcwd) then
echo "ERROR: cannot find $glmfitcwd"
exit 1;
endif
set anattype = volume;
set subject = ();
set hemi = ();
set surf = "white";
set y = ();
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 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 "--y"
set y = $glmfitcwd/$glmfitcmd[1]; shift glmfitcmd;
if(! -e $y) then
set y = $glmdir/`basename $y`
if(! -e $y) then
echo "ERROR: cannot find $y"
exit 1;
endif
endif
breaksw
case "--wls"
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
# Ignore these options that have no args
case "cmdline"
case "mri_glmfit"
case "--osgm"
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 "dods"
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"
shift glmfitcmd;
breaksw
default
echo ""
echo "WARNING: unrecognized mri_glmfit cmd option $flag"
echo ""
sleep 10
#exit 1;
breaksw
endsw
end
if($UseGRF && $anattype != volume) then
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
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
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
if($UseCache) then
set fwhm = `cat $glmdir/fwhm.dat`;
set fwhmStr = `perl -e "printf('"'%02d'"',int ( $fwhm+.5 ) ) "`
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
exit 1;
endif
echo CSD $csdCache | tee -a $LF
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
mkdir -p $glmdir/csd
@ nthJob = 1
while($nthJob <= $nJobs)
set jpad = `printf j%03d $nthJob`
echo "$nthJob/$nJobs `date`" | tee -a $LF
set cmd = (mri_glmfit --y $y $clist2 --mask $mask)
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($#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)
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)
echo "------------------------------------------------" | tee -a $LF
echo $cmd | tee -a $LF
pbsubmit $PBNodeType -c "$cmd"
sleep 5;
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)
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 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 $csdlist --mask $mask \
--cwsig $cwsig --vwsig $vwsig --sum $sum --ocn $ocn \
--oannot $oannot --annot $annot --csdpdf $csdpdf \
--cwpvalthresh $cwpvalthresh --o $msig --no-fixmni)
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;
# 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 \
--avgwf $avgwf --sum /tmp/mri_glmfit-sim.junk.$$)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
rm -f /tmp/mri_glmfit-sim.junk.$$
endif
if($anattype == volume) then
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)
if($Bonferroni) set cmd = ($cmd --bonferroni $Bonferroni)
if(! $UseGRF) set cmd = ($cmd --csdpdf $csdpdf $csdlist --vwsig $vwsig )
if($UseGRF) set cmd = ($cmd --fwhmdat $glmdir/fwhm.dat --sign $simsign
--thmin $thresh)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
# 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 \
--avgwf $avgwf --sum /tmp/mri_glmfit-sim.junk.$$)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
rm -f /tmp/mri_glmfit-sim.junk.$$
if($UseGRF) then
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
date
echo "mri_glmfit-sim done"
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 "--cache":
if( $#argv < 2) goto arg2err;
set thresh = $argv[1]; shift;
set simsign = $argv[1]; shift;
set UseCache = 1;
set DoSim = 0;
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 "--cache-dir":
if( $#argv < 1) goto arg1err;
set CacheDir = $argv[1]; shift;
breaksw
case "--cache-label":
if( $#argv < 1) goto arg1err;
set CacheLabel = $argv[1]; shift;
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 "--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;
if(`hostname` != seychelles && `hostname` != launchpad) then
hostname
echo "ERROR: must be on seychelles or launchpad to pbsubmit"
exit 1;
endif
if(`hostname` == seychelles) set PBNodeType = "-l nodes=1:opteron";
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 "--perm-force":
set PermForce = 1;
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 "--no-sim":
if($#argv < 1) goto arg1err;
set csdbase = $argv[1]; shift;
set DoSim = 0;
breaksw
case "--diag-cluster":
set DiagCluster = 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
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
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"
echo ""
echo " --glmdir glmdir"
echo " --cwp thresh : cluster p-value thresh ($cwpvalthresh)"
echo ""
echo " --cache threshold sign"
echo " --cache-dir dir : default is FREESURFER_HOME/average/mult-comp-cor"
echo " --cache-label label : default is cortex"
echo ""
echo " --grf threshold sign : -log10(p) and pos or neg (volumes only)"
echo ""
echo " --sim nulltype nsim threshold csdbase"
echo " --sim-sign sign : abs,pos,neg"
echo ""
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 " --no-sim csdbase : do not simulate, only run cluster"
echo " --seed seed : set simluation seed"
echo " --fwhm-override fwhm : override fwhm in glmdir":
echo " --uniform min max : use uniform PDF instead of gaussian":
echo ""
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 GLM-based simulations for corrections
for multiple comparisons 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 whether you want to use a pre-cached simulation or whether
you want to run a new simulation. You should use a pre-cached
simulation if: this is a surface-based analysis, you are using
fsaverage (or equivalent), you have whole-brain data, you want to
use a voxel-wise threshold of 1.3, 2.0, 2.3, 3.0, 3.3, or 4.0, and
you are going to use an mc-z simulation. Pre-cached will finish
after a few seconds, your own simulation will take hours or days.
2.a If you want to use a pre-cached, then specify "--cache threshold sign"
where threshold is the vertex-wise threshhold and sign is abs, pos,
or neg. Note: the csdbase will be set to cache.thresh.simsign (eg,
if the threshold was set to 1.3 and sign to pos, then cache.th13.pos.
Skip the remaining steps below.
2.b If you want to run your own simulation, then choose the type of
simulation (nulltype) you want to run. The choices are:
mc-z - synthesized, smoothed z-field (fast)
mc-full - fully synthesized, smoothed gaussian field (slowest)
perm - permutation on original data (probably fasted one)
3. Choose the voxel/vertex-wise threshhold used to define clusters. This
will be -log10(p), where p is the p-value. Eg, to use all voxels/vertices
with p < .01, then set the threshhold = 2.
4. Choose the number of iterations (nsim). This is usually 10,000
(which can take quite a while).
5. Select the CSD base name (csdbase). CSD stands for "Cluster
Simulation Data". There will be a text file for each contrast with
the information about the null distribution. This file will be
called csdbase-contrast.csd, so the csdbase should be something
descriptive. Eg, if you have mc-z, with threshold=2, and sign=pos,
then setting csdbase to mc-z.pos.2 would be appropriate. If the
CSD file exists, then you must delete it or run with --overwrite.
6. Select the threshold sign (--sim-sign). This is the sign of the
threshold used to create the clusters for cluster-wise correction
of multiple comparisons. Options are abs (unsigned), pos
(positive), and neg (negative). If you have an a priori hypothesis
(eg, group1 > group2), then choose positive. If you don not know
the sign, then use abs. Choosing a sign will yield stronger results
(ie, smaller p-values), but you will lose the other sign.
Now run the command:
mri_glmfit-sim --glmdir glmdir \
--cache threshold sign
or
mri_glmfit-sim --glmdir glmdir \
--sim nulltype nsim threshold csdbase \
--sim-sign sign
Example,
mri_glmfit-sim --glmdir lh.thickness.sm05 \
--cache 3 --sim-sign pos
or
mri_glmfit-sim --glmdir lh.thickness.sm05 \
--sim mc-z 10000 3 mc-z.pos.3 --sim-sign pos
This will run mri_glmfit in simulation mode. It will find all the
contrasts that you ran in the original invocation of mri_glmfit and
create CSD files for them. It will use the mask and, for mc
simulations, the FWHM found in the glmdir.
It will then run mri_surfcluster or mri_volcluster for
each contrast, using the sig volume as input and creating the
following files in each contrast directory:
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 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 stimulation but just want to run
the correction for multiple comparisons again to change the
clusterwise threshold. "csdbase" is the same that was used in the
orignial invocation to mri_glmfit-sim. Eg, if this was your original
mri_glmfit-sim command:
mri_glmfit-sim --glmdir lh.thickness.sm05 --cwp .05\
--sim mc-z 10000 3 mc-z.pos.3 --sim-sign pos
But then then you wanted to re-run with a higher CWP, then
mri_glmfit-sim --glmdir lh.thickness.sm05 --no-sim mc-z.pos.3 --cwp .1
--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.
_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
The information in this e-mail is intended only for the person to whom it is
addressed. If you believe this e-mail was sent to you in error and the e-mail
contains patient information, please contact the Partners Compliance HelpLine at
http://www.partners.org/complianceline . If the e-mail was sent to you in error
but does not contain patient information, please contact the sender and properly
dispose of the e-mail.
