User Tools

Site Tools


biac:cluster:examples:fdr

FSL False Discovery Rate Example Script

The example script below will run FSL's FDR calculation on a 3rd level cope image. The output will be a p-map (activation map of P-values) masked by the FDR threshold calculated (based on the degrees-of-freedom in the 3rd level analysis), as well as text output file that will list the various thresholds that were calculated. This script is based on information taken from the FDR User Guide.

fsl_fdr_allcopes.sh

This script will loop through all copes sending each one to a different computing node for speed.

#!/bin/bash
 
FDRRATE=0.05
 
for i in 1 2 3 4 5 6 7 8 9
do
  qsub -v EXPERIMENT=Dummy.01 qsub_fsl_fdr.sh $i $FDRRATE
done

qsub_fsl_fdr.sh

Example usage:

[user@head ~]$ qsub -v EXPERIMENT=Dummy.01 qsub_fsl_fdr.sh 1 0.05

Copy and edit script below:

#!/bin/sh
 
# This is a BIAC template script for jobs on the cluster
# You have to provide the Experiment on command line  
# when you submit the job the cluster.
#
# >  qsub -v EXPERIMENT=Dummy.01  script.sh args
#
# There are 2 USER sections 
#  1. USER DIRECTIVE: If you want mail notifications when
#     your job is completed or fails you need to set the 
#     correct email address.
#		   
#  2. USER SCRIPT: Add the user script in this section.
#     Within this section you can access your experiment 
#     folder using $EXPERIMENT. All paths are relative to this variable
#     eg: $EXPERIMENT/Data $EXPERIMENT/Analysis	
#     By default all terminal output is routed to the " Analysis "
#     folder under the Experiment directory i.e. $EXPERIMENT/Analysis
#     To change this path, set the OUTDIR variable in this section
#     to another location under your experiment folder
#     eg: OUTDIR=$EXPERIMENT/Analysis/GridOut 	
#     By default on successful completion the job will return 0
#     If you need to set another return code, set the RETURNCODE
#     variable in this section. To avoid conflict with system return 
#     codes, set a RETURNCODE higher than 100.
#     eg: RETURNCODE=110
#     Arguments to the USER SCRIPT are accessible in the usual fashion
#     eg:  $1 $2 $3
# The remaining sections are setup related and don't require
# modifications for most scripts. They are critical for access
# to your data  	 
 
# --- BEGIN GLOBAL DIRECTIVE -- 
#$ -S /bin/sh
#$ -o $HOME/$JOB_NAME.$JOB_ID.out
#$ -e $HOME/$JOB_NAME.$JOB_ID.out
#$ -m ea
# -- END GLOBAL DIRECTIVE -- 
 
# -- BEGIN PRE-USER --
#Name of experiment whose data you want to access 
EXPERIMENT=${EXPERIMENT:?"Experiment not provided"}
 
EXPERIMENT=`findexp $EXPERIMENT`
EXPERIMENT=${EXPERIMENT:?"Returned NULL Experiment"}
 
if [ $EXPERIMENT = "ERROR" ]
then
	exit 32
else 
#Timestamp
echo "----JOB [$JOB_NAME.$JOB_ID] START [`date`] on HOST [$HOSTNAME]----" 
# -- END PRE-USER --
# **********************************************************
 
# -- BEGIN USER DIRECTIVE --
# Send notifications to the following address
#$ -M yourname@somewhere.edu
 
# -- END USER DIRECTIVE --
 
# -- BEGIN USER SCRIPT --
# User script goes here
 
COPE=$1
FDRRATE=$2
 
STATSDIR=$EXPERIMENT/Analysis/FSL/3rdLevel/GroupTrialtypeA.gfeat/cope1.feat/stats
OUTDIR=$STATSDIR
 
# Convert T-values to P-values using the cope and variance cope images
cd $STATSDIR
ttologp -logpout logp${COPE} varcope${COPE} cope${COPE} `cat dof`
 
# Caclulate the p-value image
fslmaths logp${COPE} -exp p${COPE}
 
# Calculate the FDR from the p-value image
fdr -i p${COPE} -m ../mask -q $FDRRATE > fdr${COPE}
 
# Get the FDR threshold
echo "FDR (rate) is" $FDRRATE
FDRTHRESH=`awk 'NR==2' fdr${COPE}`
echo "FDR threshold is" $FDRTHRESH
FDRTHRESH=`awk -v a=$FDRTHRESH -v b="1" 'BEGIN{print (b - a)}'`
echo "1 minus FDR threshold is" $FDRTHRESH
 
# Calculate a p-value image using the FDR threshold 
fslmaths p${COPE} -mul -1 -add 1 -thr $FDRTHRESH -mas ../mask thresh_fdr_1_minus_p${COPE}
 
# Clean up unneccesary files
rm -f logp${COPE}.nii.gz p${COPE}.nii.gz fdr${COPE}
 
# -- END USER SCRIPT -- #
 
# **********************************************************
# -- BEGIN POST-USER -- 
echo "----JOB [$JOB_NAME.$JOB_ID] STOP [`date`]----" 
OUTDIR=${OUTDIR:-$EXPERIMENT/Analysis}
mv $HOME/$JOB_NAME.$JOB_ID.out $OUTDIR/$JOB_NAME.cope${COPE}.$JOB_ID.out	 
RETURNCODE=${RETURNCODE:-0}
exit $RETURNCODE
fi
# -- END POST USER--
biac/cluster/examples/fdr.txt · Last modified: 2023/02/23 18:43 (external edit)