Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#!/bin/sh
# dual_regression - take group-ICA maps (etc) and get subject-specific versions of them (and associated timecourses)
#
# Stephen Smith, Christian Beckmann, Janine Bijsterbosch, Sam Harrison, FMRIB Image Analysis Group
#
# Copyright (C) 2011-2019 University of Oxford
#
# SHCOPYRIGHT
export LC_ALL=C
Usage() {
cat <<EOF
dual_regression v0.6
***NOTE*** ORDER OF COMMAND-LINE ARGUMENTS IS DIFFERENT FROM PREVIOUS VERSION
Usage: dual_regression <group_IC_maps> <des_norm> <design.mat> <design.con> <n_perm> [--thr] <output_directory> <input1> <input2> <input3> .........
e.g. dual_regression groupICA.gica/groupmelodic.ica/melodic_IC 1 design.mat design.con 500 0 grot \`cat groupICA.gica/.filelist\`
<group_IC_maps_4D> 4D image containing spatial IC maps (melodic_IC) from the whole-group ICA analysis
<des_norm> 0 or 1 (1 is recommended). Whether to variance-normalise the timecourses used as the stage-2 regressors
<design.mat> Design matrix for final cross-subject modelling with randomise
<design.con> Design contrasts for final cross-subject modelling with randomise
<n_perm> Number of permutations for randomise; set to 1 for just raw tstat output, set to 0 to not run randomise at all.
[--thr] Perform thresholded dual regression to obtain unbiased timeseries for connectomics analyses (e.g., with FSLnets)
<output_directory> This directory will be created to hold all output and logfilesg
<input1> <input2> ... List all subjects' preprocessed, standard-space 4D datasets
<design.mat> <design.con> can be replaced with just
-1 for group-mean (one-group t-test) modelling.
If you need to add other randomise options then edit the line after "EDIT HERE" in the dual_regression script
EOF
exit 1
}
############################################################################
[ "$6" = "" ] && Usage
ORIG_COMMAND=$*
ICA_MAPS=`${FSLDIR}/bin/remove_ext $1` ; shift
DES_NORM=--des_norm
if [ $1 = 0 ] ; then
DES_NORM=""
fi ; shift
if [ $1 = "-1" ] ; then
DESIGN="-1"
shift
else
dm=$1
dc=$2
DESIGN="-d $1 -t $2"
shift 2
fi
NPERM=$1 ; shift
NAF2=0
if [ $1 = "--thr" ] ; then
NAF2=1
shift
fi
OUTPUT=`${FSLDIR}/bin/remove_ext $1` ; shift
while [ _$1 != _ ] ; do
INPUTS="$INPUTS `${FSLDIR}/bin/remove_ext $1`"
shift
done
############################################################################
mkdir $OUTPUT
LOGDIR=${OUTPUT}/scripts+logs
mkdir $LOGDIR
echo $ORIG_COMMAND > $LOGDIR/command
if [ "$DESIGN" != -1 ] ; then
/bin/cp $dm $OUTPUT/design.mat
/bin/cp $dc $OUTPUT/design.con
fi
# Key sizes
N_ICs=`$FSLDIR/bin/fslnvols $ICA_MAPS`
echo "creating common mask"
j=0
for i in $INPUTS ; do
echo "$FSLDIR/bin/fslmaths $i -Tstd -bin ${OUTPUT}/mask_`${FSLDIR}/bin/zeropad $j 5` -odt char" >> ${LOGDIR}/drA
j=`echo "$j 1 + p" | dc -`
done
ID_drA=`$FSLDIR/bin/fsl_sub -T 10 -N mask_generation1 -l $LOGDIR -t ${LOGDIR}/drA`
cat <<EOF > ${LOGDIR}/drB
#!/bin/sh
\$FSLDIR/bin/fslmerge -t ${OUTPUT}/maskALL \`\$FSLDIR/bin/imglob ${OUTPUT}/mask_*\`
\$FSLDIR/bin/fslmaths $OUTPUT/maskALL -Tmin $OUTPUT/mask
\$FSLDIR/bin/imrm $OUTPUT/mask_*
EOF
chmod a+x ${LOGDIR}/drB
ID_drB=`$FSLDIR/bin/fsl_sub -j $ID_drA -T 5 -N mask_generation2 -l $LOGDIR ${LOGDIR}/drB`
echo "doing the dual regressions"
j=0
for i in $INPUTS ; do
s=subject`${FSLDIR}/bin/zeropad $j 5`
echo "$FSLDIR/bin/fsl_glm -i $i -d $ICA_MAPS -o $OUTPUT/dr_stage1_${s}.txt --demean -m $OUTPUT/mask ; \
$FSLDIR/bin/fsl_glm -i $i -d $OUTPUT/dr_stage1_${s}.txt -o $OUTPUT/dr_stage2_$s --out_z=$OUTPUT/dr_stage2_${s}_Z --demean -m $OUTPUT/mask $DES_NORM ; \
$FSLDIR/bin/fslsplit $OUTPUT/dr_stage2_$s $OUTPUT/dr_stage2_${s}_ic" >> ${LOGDIR}/drC
j=`echo "$j 1 + p" | dc -`
done
ID_drC=`$FSLDIR/bin/fsl_sub -j $ID_drB -T 30 -N dual_regression -l $LOGDIR -t ${LOGDIR}/drC`
# For thresholded maps, try and normalise such that null is N(0,1)
# Because of the presence of signal in the tails, we use median and IQR as
# robust estimators for the mean & standard deviation
# https://en.wikipedia.org/wiki/Robust_measures_of_scale
if [ $NAF2 -eq 1 ] ; then
echo "doing thresholded dual regression"
j=0
for i in $INPUTS ; do
s=subject`${FSLDIR}/bin/zeropad $j 5`
STAGE2_MAPS=${OUTPUT}/dr_stage2_${s}
STAGE4_MAPS=${OUTPUT}/dr_stage4_${s}
echo "$FSLDIR/bin/fslsplit ${STAGE2_MAPS} ${STAGE2_MAPS}_split -t; \
for im in \$($FSLDIR/bin/imglob ${STAGE2_MAPS}_split*); do \
mean=\$($FSLDIR/bin/fslstats \${im} -p 50); \
std=\$($FSLDIR/bin/fslstats \${im} -p 25 -p 75 | awk '{ print (\$2 - \$1) / 1.349 }'); \
$FSLDIR/bin/fslmaths \${im} -sub \${mean} -div \${std} \${im}_norm; \
done; \
$FSLDIR/bin/fslmerge -t ${STAGE4_MAPS}_norm \$(imglob ${STAGE2_MAPS}_split*_norm*); \
$FSLDIR/bin/imrm \$(imglob ${STAGE2_MAPS}_split*); \
$FSLDIR/bin/fslmaths ${STAGE4_MAPS}_norm -abs -thr 2 -bin ${STAGE4_MAPS}_mask; \
$FSLDIR/bin/fslmaths ${STAGE4_MAPS}_norm -mul ${STAGE4_MAPS}_mask ${STAGE4_MAPS}_thresh; \
$FSLDIR/bin/fsl_glm -i ${i} -d ${STAGE4_MAPS}_thresh -o ${OUTPUT}/dr_stage4_${s}.txt --demean -m ${OUTPUT}/mask; \
" >> ${LOGDIR}/drD
ID_drD=`$FSLDIR/bin/fsl_sub -j $ID_drC -N thresholdedDR -l $LOGDIR -t ${LOGDIR}/drD`
fi
echo "sorting maps and running randomise"
j=0
jj=`$FSLDIR/bin/zeropad $j 4`
RAND=""
if [ $NPERM -eq 1 ] ; then
RAND="$FSLDIR/bin/randomise -i $OUTPUT/dr_stage2_ic$jj -o $OUTPUT/dr_stage3_ic$jj -m $OUTPUT/mask $DESIGN -n 1 -V -R"
fi
if [ $NPERM -gt 1 ] ; then
# EDIT HERE
RAND="$FSLDIR/bin/randomise -i $OUTPUT/dr_stage2_ic$jj -o $OUTPUT/dr_stage3_ic$jj -m $OUTPUT/mask $DESIGN -n $NPERM -T -V"
fi
echo "$FSLDIR/bin/fslmerge -t $OUTPUT/dr_stage2_ic$jj \`\$FSLDIR/bin/imglob $OUTPUT/dr_stage2_subject*_ic${jj}.*\` ; \
$FSLDIR/bin/imrm \`\$FSLDIR/bin/imglob $OUTPUT/dr_stage2_subject*_ic${jj}.*\` ; $RAND" >> ${LOGDIR}/drE
j=`echo "$j 1 + p" | dc -`
done
ID_drE=`$FSLDIR/bin/fsl_sub -j $ID_drC -T 60 -N randomise -l $LOGDIR -t ${LOGDIR}/drE`