From 6e67a828185fd00fa270fd8b71de52bcd875dbc4 Mon Sep 17 00:00:00 2001
From: Fidel Alfaro Almagro <falmagro@fmrib.ox.ac.uk>
Date: Fri, 24 Feb 2023 11:48:35 +0000
Subject: [PATCH] Improving IDPs scripts

---
 bip/data/FileTree.tree                        |   4 +-
 bip/data/tfMRI_group/groupMask1.nii.gz        |   3 +
 bip/data/tfMRI_group/groupMask2.nii.gz        |   3 +
 bip/data/tfMRI_group/groupMask5.nii.gz        |   3 +
 bip/data/tfMRI_group/groupMask5a.nii.gz       |   3 +
 bip/pipelines/IDPs_gen/IDP_SWI_T2star.py      |  42 ++++++
 bip/pipelines/IDPs_gen/IDP_T1_FIRST_vols.py   |  10 +-
 .../IDPs_gen/IDP_T1_GM_parcellation.py        |   4 +-
 bip/pipelines/IDPs_gen/IDP_T1_align_to_std.py |  54 ++++++++
 bip/pipelines/IDPs_gen/IDP_T1_noise_ratio.py  |  50 ++++++++
 bip/pipelines/IDPs_gen/IDP_all_align_to_T1.py |  55 ++++++++
 bip/pipelines/IDPs_gen/IDP_diff_autoPtx.py    | 120 ++++++++++++++++++
 .../IDPs_gen/IDP_func_task_activation.py      | 114 +++++++++++++++++
 bip/pipelines/IDPs_gen/IDPs_gen.py            |  35 +++++
 bip/pipelines/fMRI_rest/rfMRI_melodic.py      |  23 ++--
 bip/pipelines/fMRI_task/tfMRI_feat.py         |  34 ++++-
 bip/utils/log_utils.py                        |  10 +-
 17 files changed, 540 insertions(+), 27 deletions(-)
 create mode 100644 bip/data/tfMRI_group/groupMask1.nii.gz
 create mode 100644 bip/data/tfMRI_group/groupMask2.nii.gz
 create mode 100644 bip/data/tfMRI_group/groupMask5.nii.gz
 create mode 100644 bip/data/tfMRI_group/groupMask5a.nii.gz
 create mode 100755 bip/pipelines/IDPs_gen/IDP_SWI_T2star.py
 create mode 100755 bip/pipelines/IDPs_gen/IDP_T1_align_to_std.py
 create mode 100755 bip/pipelines/IDPs_gen/IDP_T1_noise_ratio.py
 create mode 100755 bip/pipelines/IDPs_gen/IDP_all_align_to_T1.py
 create mode 100755 bip/pipelines/IDPs_gen/IDP_diff_autoPtx.py
 create mode 100755 bip/pipelines/IDPs_gen/IDP_func_task_activation.py

diff --git a/bip/data/FileTree.tree b/bip/data/FileTree.tree
index 7b6cffb..d7f4238 100644
--- a/bip/data/FileTree.tree
+++ b/bip/data/FileTree.tree
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:0d593e72e427400f23ec2829c5353cea47df1516fd3908e3a1b6c0c94f5692e6
-size 22258
+oid sha256:30abf95478fe8611d0162a4bce17d4cabae1471ea0d322398d1071a8417e0009
+size 23361
diff --git a/bip/data/tfMRI_group/groupMask1.nii.gz b/bip/data/tfMRI_group/groupMask1.nii.gz
new file mode 100644
index 0000000..09a1f59
--- /dev/null
+++ b/bip/data/tfMRI_group/groupMask1.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ae191189dfb9dd61c615695524b609712b891091e0bb9268627b14334bf09191
+size 12596
diff --git a/bip/data/tfMRI_group/groupMask2.nii.gz b/bip/data/tfMRI_group/groupMask2.nii.gz
new file mode 100644
index 0000000..763c578
--- /dev/null
+++ b/bip/data/tfMRI_group/groupMask2.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:67aad7cba913142158726cc90fcffda64bf5b4d50f0e87e3362a1f537ae9d6db
+size 12315
diff --git a/bip/data/tfMRI_group/groupMask5.nii.gz b/bip/data/tfMRI_group/groupMask5.nii.gz
new file mode 100644
index 0000000..dabbafd
--- /dev/null
+++ b/bip/data/tfMRI_group/groupMask5.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bad580a85019111b481be06bde3704ee8866d75dabe89ff390a65cae95ba0532
+size 5642
diff --git a/bip/data/tfMRI_group/groupMask5a.nii.gz b/bip/data/tfMRI_group/groupMask5a.nii.gz
new file mode 100644
index 0000000..610f4b8
--- /dev/null
+++ b/bip/data/tfMRI_group/groupMask5a.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:eea197dc41c09f8f907b0743a6a2c13dd0e72f40d0e3a4d4a3f4ceb7ba23ba07
+size 4082
diff --git a/bip/pipelines/IDPs_gen/IDP_SWI_T2star.py b/bip/pipelines/IDPs_gen/IDP_SWI_T2star.py
new file mode 100755
index 0000000..52d66e1
--- /dev/null
+++ b/bip/pipelines/IDPs_gen/IDP_SWI_T2star.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+#
+# IDP_SWI_T2star.py - Generating IDP file with metrics of T2star
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import json
+import logging
+from fsl import wrappers
+from pipe_tree import In, Out, Ref
+from bip.utils.log_utils import redirect_logging
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T2star_to_T1:               In,
+        T1_first_all_fast_firstseg: In(optional=True),        
+        logs_dir:                   Ref,
+        IDP_SWI_T2star:             Out):
+    
+    with redirect_logging('IDP_SWI_T2star', outdir=logs_dir):
+
+        result = ("NaN " * 14).strip()
+
+        if os.path.exists(T1_first_all_fast_firstseg):
+            v=wrappers.fslstats(T2star_to_T1, K=T1_first_all_fast_firstseg).p(50).run()
+
+            # Check that the outputs are OK
+            if len(v) == 58:
+                # indices that we are interested in
+                ind = [9, 48, 10, 49, 11, 50, 12, 51, 16, 52, 17, 53, 25, 57]
+                result = [str(int(v[x])) for x in ind]
+                result = " ".join(result)
+
+        with open(IDP_SWI_T2star, 'wt', encoding="utf-8") as f:
+            f.write(f'{result}\n')
\ No newline at end of file
diff --git a/bip/pipelines/IDPs_gen/IDP_T1_FIRST_vols.py b/bip/pipelines/IDPs_gen/IDP_T1_FIRST_vols.py
index 5c5b166..d261eb0 100755
--- a/bip/pipelines/IDPs_gen/IDP_T1_FIRST_vols.py
+++ b/bip/pipelines/IDPs_gen/IDP_T1_FIRST_vols.py
@@ -28,10 +28,12 @@ def run(ctx,
 
         if os.path.exists(T1_first_all_fast_firstseg):
             v=wrappers.fslstats(T1_first_all_fast_firstseg).H(58,0.5,58.5).run()
-            # indices that we are interested in
-            ind = [9, 48, 10, 49, 11, 50, 12, 51, 16, 52, 17, 53, 25, 57, 15]
-            result = [str(int(v[x])) for x in ind]
-            result = " ".join(result)
+            # Check that the outputs are OK
+            if len(v) == 58:
+                # indices that we are interested in
+                ind = [9, 48, 10, 49, 11, 50, 12, 51, 16, 52, 17, 53, 25, 57,15]
+                result = [str(int(v[x])) for x in ind]
+                result = " ".join(result)
 
 
         with open(IDP_T1_FIRST_vols, 'wt', encoding="utf-8") as f:
diff --git a/bip/pipelines/IDPs_gen/IDP_T1_GM_parcellation.py b/bip/pipelines/IDPs_gen/IDP_T1_GM_parcellation.py
index 2fe4efd..e7a9fbe 100755
--- a/bip/pipelines/IDPs_gen/IDP_T1_GM_parcellation.py
+++ b/bip/pipelines/IDPs_gen/IDP_T1_GM_parcellation.py
@@ -22,6 +22,7 @@ def run(ctx,
         T1:                      In,
         T1_to_MNI_warp_coef_inv: In,
         T1_fast_pve_1:           In,
+        tmp_dir:                 Ref,
         logs_dir:                Ref,
         IDP_T1_GM_parcellation:  Out):
     
@@ -35,8 +36,7 @@ def run(ctx,
         wrappers.applywarp(src=GMatlas, ref=T1, w=T1_to_MNI_warp_coef_inv,
                            out=GMatlas_to_T1, interp='nn')
 
-        wrappers.fslstats(fMRI_SNR).l(0.1).p(50).run()
-        vals = wrappers.fslstats(T1_brain_pve_1, K=GMatlas_to_T1).m.v.run()
+        vals = wrappers.fslstats(T1_fast_pve_1, K=GMatlas_to_T1).m.v.run()
         result = " ".join([str(x) for x in vals[:,0] * vals[:,1]])
 
         with open(IDP_T1_GM_parcellation, 'wt', encoding="utf-8") as f:
diff --git a/bip/pipelines/IDPs_gen/IDP_T1_align_to_std.py b/bip/pipelines/IDPs_gen/IDP_T1_align_to_std.py
new file mode 100755
index 0000000..f6de927
--- /dev/null
+++ b/bip/pipelines/IDPs_gen/IDP_T1_align_to_std.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+#
+# IDP_T1_align_to_std.py - Generating IDP file with alignment to std metrics.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import logging
+from fsl import wrappers
+from pipe_tree import In, Out, Ref
+from bip.utils.log_utils import redirect_logging, tempdir
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_brain:             In,
+        T1_to_MNI_linear_mat: In,
+        T1_brain_to_MNI:      In,
+        T1_to_MNI_warp_jac:   In,
+        tmp_dir:              Ref,
+        logs_dir:             Ref,
+        IDP_T1_align_to_std:  Out):
+    
+    with redirect_logging('IDP_T1_align_to_std', outdir=logs_dir),\
+        tempdir(tmp_dir):
+
+        tmp_jac = tmp_dir + '/tmpjac.nii.gz'
+        tmp_mat = tmp_dir + '/tmp_mat.mat' 
+
+        MC=ctx.FSLDIR + '/etc/flirtsch/measurecost1.sch'
+        MNI152_T1_1mm_brain = ctx.get_standard("MNI152_T1_1mm_brain.nii.gz")
+        MNI152_T1_1mm_brain_mask = ctx.get_standard("MNI152_T1_1mm_brain_mask.nii.gz")
+
+        costs1 = wrappers.flirt(src=T1_brain, ref=MNI152_T1_1mm_brain, 
+                                refweight=MNI152_T1_1mm_brain_mask, 
+                                init=T1_to_MNI_linear_mat, 
+                                schedule=MC, omat=tmp_mat).stdout[0].strip()
+        val1 = costs1.split()[0]
+
+        costs2 = wrappers.flirt(src=T1_brain_to_MNI, ref=MNI152_T1_1mm_brain, 
+                                refweight=MNI152_T1_1mm_brain_mask, 
+                                schedule=MC, omat=tmp_mat).stdout[0].strip()
+        val2 = costs2.split()[0]
+
+        wrappers.fslmaths(T1_to_MNI_warp_jac).sub(1).sqr().run(tmp_jac)
+        val3 = wrappers.fslstats(tmp_jac, K=MNI152_T1_1mm_brain_mask).m.run()
+
+        with open(IDP_T1_align_to_std, 'wt', encoding="utf-8") as f:
+            f.write(f'{val1} {val2} {val3}\n')
diff --git a/bip/pipelines/IDPs_gen/IDP_T1_noise_ratio.py b/bip/pipelines/IDPs_gen/IDP_T1_noise_ratio.py
new file mode 100755
index 0000000..82a2921
--- /dev/null
+++ b/bip/pipelines/IDPs_gen/IDP_T1_noise_ratio.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+#
+# IDP_T1_noise_ratio.py - Generating IDP file with measure of noise ratio of T1.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import json
+import logging
+from fsl import wrappers
+from pipe_tree import In, Out, Ref
+from bip.utils.log_utils import redirect_logging, tempdir
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_fast_pveseg:    In,
+        T1:                 In,
+        tmp_dir:            Ref,
+        logs_dir:           Ref,
+        IDP_T1_noise_ratio: Out):
+    
+    with redirect_logging('IDP_T1_noise_ratio', outdir=logs_dir),\
+        tempdir(tmp_dir):
+
+        tmp_SNR = tmp_dir + '/tmp_SNR.nii.gz'
+
+        wrappers.fslmaths(T1_fast_pveseg).thr(2).uthr(2).ero().run(tmp_SNR)
+        grey  = wrappers.fslstats(T1).k(tmp_SNR).m.run()
+        wrappers.fslmaths(T1_fast_pveseg).thr(3).uthr(3).ero().run(tmp_SNR)
+        white = wrappers.fslstats(T1).k(tmp_SNR).m.run()
+
+        #TODO: The roundings are only there to mimic the behaviour of
+        #      previous version. We do not have to do it in next version
+        brain    = int(round((grey + white) / 2))
+        contrast = int(round(white - grey))
+        thresh   = int(round(brain / 10))
+
+        noise = wrappers.fslstats(T1).l(0.001).u(thresh).s.run()
+
+        SNR = noise / brain
+        CNR = noise / contrast
+
+        with open(IDP_T1_noise_ratio, 'wt', encoding="utf-8") as f:
+            f.write(f'{SNR} {CNR}\n')
diff --git a/bip/pipelines/IDPs_gen/IDP_all_align_to_T1.py b/bip/pipelines/IDPs_gen/IDP_all_align_to_T1.py
new file mode 100755
index 0000000..f6155d4
--- /dev/null
+++ b/bip/pipelines/IDPs_gen/IDP_all_align_to_T1.py
@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+#
+# IDP_all_align_to_T1.py - Generating IDP file with alignment to T1 metric.
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import logging
+from fsl import wrappers
+from pipe_tree import In, Out, Ref
+from bip.utils.log_utils import redirect_logging, tempdir
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_brain:                   In,
+        T1_brain_mask:              In,
+        T2_FLAIR_brain:             In(optional=True),
+        fieldmap_iout_to_T1:        In(optional=True),
+        SWI_to_T1:                  In(optional=True),
+        rfMRI_example_func2highres: In(optional=True),
+        tfMRI_example_func2highres: In(optional=True),
+        tmp_dir:                    Ref,
+        logs_dir:                   Ref,
+        IDP_all_align_to_T1:        Out):
+
+    with redirect_logging('IDP_all_align_to_T1', outdir=logs_dir),\
+        tempdir(tmp_dir):
+
+        tmp_mat = tmp_dir + '/tmp_mat.mat' 
+
+        result=""
+
+        MC=ctx.FSLDIR + '/etc/flirtsch/measurecost1.sch'
+
+        for file_name in [T2_FLAIR_brain, fieldmap_iout_to_T1, SWI_to_T1, 
+                         rfMRI_example_func2highres,tfMRI_example_func2highres]:
+
+            if os.path.exists(file_name):
+                costs1 = wrappers.flirt(src=file_name, ref=T1_brain, 
+                                        refweight=T1_brain_mask, schedule=MC, 
+                                        omat=tmp_mat).stdout[0].strip()
+                result += " " + str(costs1.split()[0])
+            else:
+                result += " NaN"
+
+        result = result.strip()
+
+        with open(IDP_all_align_to_T1, 'wt', encoding="utf-8") as f:
+            f.write(f'{result}\n')
diff --git a/bip/pipelines/IDPs_gen/IDP_diff_autoPtx.py b/bip/pipelines/IDPs_gen/IDP_diff_autoPtx.py
new file mode 100755
index 0000000..786285e
--- /dev/null
+++ b/bip/pipelines/IDPs_gen/IDP_diff_autoPtx.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+#
+# IDP_diff_autoPtx.py - Generating IDP file with autoPtx metrics 
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import json
+import logging
+from pipe_tree import In, Out, Ref
+from bip.utils.log_utils import redirect_logging
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_orig:            In,
+        logs_dir:           Ref,
+        IDP_diff_autoPtx: Out):
+    
+    with redirect_logging('IDP_diff_autoPtx', outdir=logs_dir):
+
+
+set -x
+
+origDir=`pwd`
+scriptName=`basename "$0"`
+direc=$1
+
+cd $direc
+
+basedMRI="dMRI"
+#if [ -d $basedMRI/unusable ] ; then
+#    basedMRI="$basedMRI/unusable"
+#elif [ -d $basedMRI/incompatible ] ; then
+#    basedMRI="$basedMRI/incompatible"
+#fi
+
+#Setting the string of NaN in case there is a problem.
+numVars="27"
+nanResult="";
+for i in $(seq 1 $numVars) ; do 
+    nanResult="NaN $nanResult" ; 
+done
+
+result="" 
+
+allgood=1
+for i in FA MD MO L1 L2 L3 ICVF OD ISOVF ; do
+  if [ ! -f $basedMRI/autoptx_preproc/APTX_${i}.txt ] || [ `cat $basedMRI/autoptx_preproc/APTX_${i}.txt | wc -w` != 27 ] ; then
+    allgood=0
+  fi
+done
+
+if [ $allgood = 0 ] ; then
+  if [ -f $basedMRI/autoptx_preproc/tracts/mcp/tracts/tractsNorm.nii.gz ] ; then
+
+    ntracts=`echo $basedMRI/autoptx_preproc/tracts/*/tracts/tractsNorm.nii.gz | wc -w`
+    if [ $ntracts != 27 ] ; then
+      >&2 echo only found $ntracts tracts
+    fi
+
+    # hack to ameliorate partial misalignment of mcp
+    fslmaths $FSLDIR/data/standard/FMRIB58_FA_1mm -thr 1000 -bin -ero -ero -mul $basedMRI/autoptx_preproc/tracts/mcp/tracts/tractsNorm $basedMRI/autoptx_preproc/tracts/mcp/tracts/tractsNorm -odt float
+
+    for p in ar_l ar_r atr_l atr_r cgc_l cgc_r cgh_l cgh_r cst_l cst_r fma fmi ifo_l ifo_r ilf_l ilf_r mcp ml_l ml_r ptr_l ptr_r slf_l slf_r str_l str_r unc_l unc_r ; do
+      if [ -f $basedMRI/autoptx_preproc/tracts/$p/tracts/tractsNorm.nii.gz ] ; then
+        t=`fslstats $basedMRI/autoptx_preproc/tracts/$p/tracts/tractsNorm -P 99`
+        fslmaths $basedMRI/autoptx_preproc/tracts/$p/tracts/tractsNorm -min $t -mul 1e6 -div $t -thr 1e5 $basedMRI/autoptx_preproc/tmp_$p
+      fi
+    done
+
+    cad=""
+
+    fslmerge -t $basedMRI/autoptx_preproc/tmp_all \
+         $basedMRI/autoptx_preproc/tmp_{ar_l,ar_r,atr_l,atr_r,cgc_l,cgc_r,cgh_l,cgh_r,cst_l,cst_r,fma,fmi,ifo_l,ifo_r,ilf_l,ilf_r,mcp,ml_l,ml_r,ptr_l,ptr_r,slf_l,slf_r,str_l,str_r,unc_l,unc_r}.nii.gz
+    tractnorm=`fslstats -t $basedMRI/autoptx_preproc/tmp_all -m`
+
+    if [ -f $basedMRI/autoptx_preproc/tmp_all.nii.gz ] ; then
+        for i in FA MD MO L1 L2 L3 ICVF OD ISOVF ; do
+          /bin/rm -f $basedMRI/autoptx_preproc/APTX_${i}.txt
+          if [ -f $basedMRI/TBSS/stats/all_${i}.nii.gz ] ; then
+              fslmaths $basedMRI/autoptx_preproc/tmp_all -mul $basedMRI/TBSS/stats/all_${i} $basedMRI/autoptx_preproc/tmp
+              tractmeas=`fslstats -t $basedMRI/autoptx_preproc/tmp -m | sed 's/-/_/g'`
+              jmax=`fslnvols $basedMRI/autoptx_preproc/tmp_all`
+              for (( j=1; j <= $jmax; ++j )) ; do
+                tm=`echo $tractmeas | cut -d " " -f$j`
+                tn=`echo $tractnorm | cut -d " " -f$j`
+                echo "10 k $tm $tn / p" | dc - >> $basedMRI/autoptx_preproc/APTX_${i}.txt
+              done
+          fi
+        done
+        rm $basedMRI/autoptx_preproc/tmp_*.nii.gz
+    fi
+  else
+    >&2 echo "mcp tract missing"
+  fi
+fi
+
+for i in FA MD MO L1 L2 L3 ICVF OD ISOVF ; do
+  if [ -f $basedMRI/autoptx_preproc/APTX_${i}.txt ] && [ `cat $basedMRI/autoptx_preproc/APTX_${i}.txt | wc -w` = 27 ] ; then
+    miniResult=`cat $basedMRI/autoptx_preproc/APTX_${i}.txt`
+  else
+    miniResult="$nanResult"
+  fi
+  result="$result $miniResult"
+done
+
+mkdir -p IDP_files
+
+echo $result > IDP_files/$scriptName.txt
+echo $result
+
+cd $origDir
+
+
diff --git a/bip/pipelines/IDPs_gen/IDP_func_task_activation.py b/bip/pipelines/IDPs_gen/IDP_func_task_activation.py
new file mode 100755
index 0000000..084d16c
--- /dev/null
+++ b/bip/pipelines/IDPs_gen/IDP_func_task_activation.py
@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+#
+# IDP_func_task_activation.py - Generating IDP file with tfMRI activation
+#
+# Author: Fidel Alfaro Almagro <fidel.alfaroalmagro@ndcn.ox.ac.uk>
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+# Author: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
+#
+# pylint: disable=C0103,E0602,C0114,C0115,C0116,R0913,R0914,R0915
+#
+
+import os
+import json
+import shutil
+import logging
+from fsl import wrappers
+from pipe_tree import In, Out, Ref
+from bip.utils.log_utils import redirect_logging
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_to_MNI_warp:            In,
+        T1_to_MNI_warp_coef_inv:   In,
+        tfMRI_feat:                In,
+        tfMRI_cope1:               In,
+        tfMRI_cope2:               In,
+        tfMRI_cope5:               In,
+        tfMRI_zstat1:              In,
+        tfMRI_zstat2:              In,
+        tfMRI_zstat5:              In,
+        highres2standard_warp:     Ref,
+        highres2standard_warp_inv: Ref,
+        logs_dir:                  Ref,
+        tfMRI_featquery_1_dir:     Out,
+        tfMRI_featquery_2_dir:     Out,
+        tfMRI_featquery_5_dir:     Out,
+        tfMRI_featquery_5a_dir:    Out,
+        tfMRI_featquery_1_report:  Out,
+        tfMRI_featquery_2_report:  Out,
+        tfMRI_featquery_5_report:  Out,
+        tfMRI_featquery_5a_report: Out,
+        IDP_func_task_activation:  Out):
+    
+    with redirect_logging('IDP_func_task_activation', outdir=logs_dir):
+        if not os.path.exists(highres2standard_warp):
+            os.symlink(src=os.getcwd() + "/" + T1_to_MNI_warp, 
+                       dst=highres2standard_warp)
+
+        if not os.path.exists(highres2standard_warp_inv):
+            os.symlink(src=os.getcwd() + "/" + T1_to_MNI_warp_coef_inv,
+                       dst=highres2standard_warp_inv)
+
+        for dir_name in [tfMRI_featquery_1_dir, tfMRI_featquery_2_dir,
+                         tfMRI_featquery_5_dir, tfMRI_featquery_5a_dir]:
+            if os.path.exists(dir_name):
+                shutil.rmtree(dir_name)
+
+        group_mask_1  = ctx.get_data("tfMRI_group/groupMask1.nii.gz")
+        group_mask_2  = ctx.get_data("tfMRI_group/groupMask2.nii.gz")
+        group_mask_5  = ctx.get_data("tfMRI_group/groupMask5.nii.gz")
+        group_mask_5a = ctx.get_data("tfMRI_group/groupMask5a.nii.gz")
+
+        tfMRI_feat += "/"
+
+        N_tfMRI_cope1 = tfMRI_cope1.replace(tfMRI_feat, "")
+        N_tfMRI_cope2 = tfMRI_cope2.replace(tfMRI_feat, "")
+        N_tfMRI_cope5 = tfMRI_cope5.replace(tfMRI_feat, "")
+
+        N_tfMRI_zstat1 = tfMRI_zstat1.replace(tfMRI_feat, "")
+        N_tfMRI_zstat2 = tfMRI_zstat2.replace(tfMRI_feat, "")
+        N_tfMRI_zstat5 = tfMRI_zstat5.replace(tfMRI_feat, "")
+
+        N_tfMRI_featquery_1_dir  = tfMRI_featquery_1_dir.replace(tfMRI_feat, "")
+        N_tfMRI_featquery_2_dir  = tfMRI_featquery_2_dir.replace(tfMRI_feat, "")
+        N_tfMRI_featquery_5_dir  = tfMRI_featquery_5_dir.replace(tfMRI_feat, "")
+        N_tfMRI_featquery_5a_dir = tfMRI_featquery_5a_dir.replace(tfMRI_feat,"")
+
+
+        wrappers.featquery(N_featdirs="1", featdir1=tfMRI_feat, N_stats="2", 
+                           stats1=[N_tfMRI_cope1, N_tfMRI_zstat1], 
+                           outputRootName=N_tfMRI_featquery_1_dir, 
+                           mask=group_mask_1)
+        wrappers.featquery(N_featdirs="1", featdir1=tfMRI_feat, N_stats="2", 
+                           stats1=[N_tfMRI_cope2, N_tfMRI_zstat2], 
+                           outputRootName=N_tfMRI_featquery_2_dir, 
+                           mask=group_mask_2)
+        wrappers.featquery(N_featdirs="1", featdir1=tfMRI_feat, N_stats="2", 
+                           stats1=[N_tfMRI_cope5, N_tfMRI_zstat5], 
+                           outputRootName=N_tfMRI_featquery_5_dir, 
+                           mask=group_mask_5)
+        wrappers.featquery(N_featdirs="1", featdir1=tfMRI_feat, N_stats="2", 
+                           stats1=[N_tfMRI_cope5, N_tfMRI_zstat5], 
+                           outputRootName=N_tfMRI_featquery_5a_dir, 
+                           mask=group_mask_5a)
+
+        result="" 
+        
+        for file_name in [tfMRI_featquery_1_report, tfMRI_featquery_2_report,
+                          tfMRI_featquery_5_report, tfMRI_featquery_5a_report]:       
+            with open(file_name, "r", encoding="utf-8") as f:
+                report = f.readlines()
+
+            a = round(float(report[0].split()[6]) / 100, 6)
+            b = round(float(report[0].split()[7]) / 100, 6)
+            c = round(float(report[-1].split()[6]), 6)
+            d = round(float(report[-1].split()[7]), 6)
+
+            result += " " + str(a) + " " + str(b) + " " + str(c) + " " + str(d)
+
+        result = result.strip()
+
+        with open(IDP_func_task_activation, 'wt', encoding="utf-8") as f:
+            f.write(f'{result}\n')
diff --git a/bip/pipelines/IDPs_gen/IDPs_gen.py b/bip/pipelines/IDPs_gen/IDPs_gen.py
index 89f15a6..769243a 100755
--- a/bip/pipelines/IDPs_gen/IDPs_gen.py
+++ b/bip/pipelines/IDPs_gen/IDPs_gen.py
@@ -13,9 +13,14 @@ import logging
 from bip.utils.log_utils    import redirect_logging
 from bip.pipelines.IDPs_gen import IDP_subject_ID, IDP_subject_centre
 from bip.pipelines.IDPs_gen import IDP_subject_COG_table
+from bip.pipelines.IDPs_gen import IDP_all_align_to_T1
 from bip.pipelines.IDPs_gen import IDP_T1_FIRST_vols, IDP_T1_SIENAX
+from bip.pipelines.IDPs_gen import IDP_T1_GM_parcellation, IDP_T1_SIENAX
+from bip.pipelines.IDPs_gen import IDP_T1_align_to_std, IDP_T1_noise_ratio
 from bip.pipelines.IDPs_gen import IDP_T2_FLAIR_WMH
+from bip.pipelines.IDPs_gen import IDP_SWI_T2star
 from bip.pipelines.IDPs_gen import IDP_func_head_motion, IDP_func_TSNR
+from bip.pipelines.IDPs_gen import IDP_func_task_activation
 from bip.pipelines.IDPs_gen import IDP_diff_eddy_outliers, IDP_diff_TBSS
 
 log = logging.getLogger(__name__)
@@ -32,6 +37,21 @@ def add_to_pipeline(ctx, pipe, tree, targets):
              kwargs={'ctx' : ctx})
         targets.append('IDP_subject_ID')
 
+        pipe(IDP_T1_align_to_std.run,
+             submit=dict(jobtime=200, name="BIP_IDP_T1_align_to_std_" + subj),
+             kwargs={'ctx' : ctx})
+        targets.append('IDP_T1_align_to_std')
+
+        pipe(IDP_T1_noise_ratio.run,
+             submit=dict(jobtime=200, name="BIP_IDP_T1_noise_ratio_" + subj),
+             kwargs={'ctx' : ctx})
+        targets.append('IDP_T1_noise_ratio')
+
+        pipe(IDP_all_align_to_T1.run,
+             submit=dict(jobtime=200, name="BIP_IDP_all_align_to_T1_" + subj),
+             kwargs={'ctx' : ctx})
+        targets.append('IDP_all_align_to_T1')
+
         pipe(IDP_func_head_motion.run,
              submit=dict(jobtime=200, name="BIP_IDP_func_head_motion_" + subj),
              kwargs={'ctx' : ctx})
@@ -57,11 +77,26 @@ def add_to_pipeline(ctx, pipe, tree, targets):
              kwargs={'ctx' : ctx})
         targets.append('IDP_T1_FIRST_vols')
 
+        pipe(IDP_T1_GM_parcellation.run,
+             submit=dict(jobtime=200, name="BIP_IDP_T1_GM_parcellation_" +subj),
+             kwargs={'ctx' : ctx})
+        targets.append('IDP_T1_GM_parcellation')
+
         pipe(IDP_T2_FLAIR_WMH.run,
              submit=dict(jobtime=200, name="BIP_IDP_T2_FLAIR_WMH_" + subj),
              kwargs={'ctx' : ctx})
         targets.append('IDP_T2_FLAIR_WMH')
 
+        pipe(IDP_SWI_T2star.run,
+             submit=dict(jobtime=200, name="BIP_IDP_SWI_T2star_" + subj),
+             kwargs={'ctx' : ctx})
+        targets.append('IDP_SWI_T2star')
+
+        pipe(IDP_func_task_activation.run,
+             submit=dict(jobtime=200,name="BIP_IDP_func_task_activation_"+subj),
+             kwargs={'ctx' : ctx})
+        targets.append('IDP_func_task_activation')
+
         pipe(IDP_diff_TBSS.run,
              submit=dict(jobtime=200, name="BIP_IDP_diff_TBSS_"+ subj),
              kwargs={'ctx' : ctx})
diff --git a/bip/pipelines/fMRI_rest/rfMRI_melodic.py b/bip/pipelines/fMRI_rest/rfMRI_melodic.py
index 9c048f9..8469ae8 100755
--- a/bip/pipelines/fMRI_rest/rfMRI_melodic.py
+++ b/bip/pipelines/fMRI_rest/rfMRI_melodic.py
@@ -24,6 +24,7 @@ def run(ctx,
         rfMRI_fsf:                        In,
         logs_dir:                         Ref,
         rfMRI_ica:                        Out,
+        rfMRI_example_func2highres:       Out,
         rfMRI_example_func2standard_warp: Out,
         rfMRI_example_func:               Out,
         rfMRI_standard:                   Out,
@@ -35,15 +36,15 @@ def run(ctx,
     with redirect_logging('rfMRI_melodic', outdir=logs_dir):
         wrappers.feat(fsf=rfMRI_fsf)
 
-    # The ending of this script is needed because the overwrite output
-    # function in MELODIC is not working, and therefore, the output is
-    # generated in a folder called "rfMI+.ica". Once this function works,
-    # we can remove these last 8 lines
-    name_dir_alt = rfMRI_ica.split(".")
-    name_dir_alt[-2] = name_dir_alt[-2] + "+"
-    name_dir_alt = ".".join(name_dir_alt)
+        # The ending of this script is needed because the overwrite output
+        # function in MELODIC is not working, and therefore, the output is
+        # generated in a folder called "rfMI+.ica". Once this function works,
+        # we can remove these last 8 lines
+        name_dir_alt = rfMRI_ica.split(".")
+        name_dir_alt[-2] = name_dir_alt[-2] + "+"
+        name_dir_alt = ".".join(name_dir_alt)
 
-    if os.path.exists(name_dir_alt):
-        if os.path.exists(rfMRI_ica):
-            shutil.rmtree(rfMRI_ica)
-        os.rename(name_dir_alt, rfMRI_ica)
+        if os.path.exists(name_dir_alt):
+            if os.path.exists(rfMRI_ica):
+                shutil.rmtree(rfMRI_ica)
+            os.rename(name_dir_alt, rfMRI_ica)
diff --git a/bip/pipelines/fMRI_task/tfMRI_feat.py b/bip/pipelines/fMRI_task/tfMRI_feat.py
index 1a9a26c..6b27552 100755
--- a/bip/pipelines/fMRI_task/tfMRI_feat.py
+++ b/bip/pipelines/fMRI_task/tfMRI_feat.py
@@ -18,10 +18,36 @@ from pipe_tree import In, Out, Ref
 log = logging.getLogger(__name__)
 
 def run(ctx,
-        tfMRI_fsf:                In,
-        logs_dir:                 Ref,
-        tfMRI_feat:               Out,
-        tfMRI_filtered_func_data: Out):
+        tfMRI_fsf:                  In,
+        logs_dir:                   Ref,
+        tfMRI_feat:                 Out,
+        tfMRI_mc_rel_mean:          Out,
+        tfMRI_filtered_func_data:   Out,
+        tfMRI_example_func2highres: Out,
+        tfMRI_cope1:                Out,
+        tfMRI_cope2:                Out,
+        tfMRI_cope3:                Out,
+        tfMRI_cope4:                Out,
+        tfMRI_cope5:                Out,
+        tfMRI_zstat1:               Out,
+        tfMRI_zstat2:               Out,
+        tfMRI_zstat3:               Out,
+        tfMRI_zstat4:               Out,
+        tfMRI_zstat5:               Out):
 
     with redirect_logging('tfMRI_feat', outdir=logs_dir):
+
         wrappers.feat(fsf=tfMRI_fsf)
+
+        # The ending of this script is needed because the overwrite output
+        # function in MELODIC is not working, and therefore, the output is
+        # generated in a folder called "rfMI+.ica". Once this function works,
+        # we can remove these last 8 lines
+        name_dir_alt = tfMRI_feat.split(".")
+        name_dir_alt[-2] = name_dir_alt[-2] + "+"
+        name_dir_alt = ".".join(name_dir_alt)
+
+        if os.path.exists(name_dir_alt):
+            if os.path.exists(tfMRI_feat):
+                shutil.rmtree(tfMRI_feat)
+            os.rename(name_dir_alt, tfMRI_feat)
diff --git a/bip/utils/log_utils.py b/bip/utils/log_utils.py
index ce1653c..2666ef0 100755
--- a/bip/utils/log_utils.py
+++ b/bip/utils/log_utils.py
@@ -84,7 +84,7 @@ def log_bip(logobj, level, msg):
 
 
 @contextlib.contextmanager
-def redirect_logging(name, outdir='.', level=logging.INFO):
+def redirect_logging(name, outdir='.', level=logging.INFO, verbose=True):
     """Context manager which temporarily redirects logging. The following log
     files are created:
 
@@ -120,15 +120,17 @@ def redirect_logging(name, outdir='.', level=logging.INFO):
                                          'stderr' : stdout,
                                          'cmd'    : logcmd}):
             try:
-                log.info(f'Started {name}')
-                log.info(f'Host {socket.gethostname()}')
+                if verbose:
+                    log.info(f'Started {name}')
+                    log.info(f'Host {socket.gethostname()}')
                 yield
             except Exception as e:
                 traceback.print_exception(type(e), e, e.__traceback__,
                                           file=stdout)
                 raise e
             finally:
-                log.info(f'Finished {name}')
+                if verbose:
+                    log.info(f'Finished {name}')
                 log.removeHandler(handler)
                 # Restore main log handler
                 if getattr(log, 'handler', None):
-- 
GitLab