Draft for FS sub-pipeline

bip/ext_wrappers/__init__.py

@@ -11,4 +11,6 @@
 them to be called from Python.
 """
 
-from bip.ext_wrappers.freesurfer   import  (recon_all)
+from bip.ext_wrappers.freesurfer import (recon_all,
+                                         asegstats2table,
+                                         aparcstats2table)

bip/ext_wrappers/freesurfer.py

@@ -8,24 +8,49 @@
 
 """
 
-
+import os
 import fsl.utils.assertions as asrt
-from fsl.utils.deprecated import deprecated
-from . import wrapperutils  as wutils
+from fsl.wrappers import wrapperutils as wutils
+
+
+import sys
+
+def eprint(*args, **kwargs):
+    print(*args, file=sys.stderr, **kwargs)
 
 
+@wutils.fileOrImage('infile')
 @wutils.cmdwrapper
-def recon_all(subjid, directive, **kwargs):
-    """Wrapper for the ``xfibres_gpu`` command.
+def recon_all(subjects_dir, directive, subjid, infile, FLAIR=None, **kwargs):
+    """Wrapper for the ``recon-all`` command.
+    """
+    #asrt.assertIsNifti(*infile)
 
-"""
+    os.environ["SUBJECTS_DIR"] = subjects_dir
+
+    cmd = ['recon-all', "-"+directive + " -s " + subjid ]
+    cmd += [ " -i " + subjects_dir + "/" + infile ]
 
-    valmap = {
-        'FLAIRpial' : wutils.SHOW_IF_TRUE,
-    }
+    if FLAIR:
+        cmd += [ " -FLAIR " + subjects_dir + "/" + FLAIR + " -FLAIRpial" ]
 
-    cmd = ['recon-all', "-"+directive] 
+    return cmd
+
+@wutils.cmdwrapper
+def asegstats2table(subjects_dir, mask, tablefile, subjects, **kwargs):
+    """Wrapper for the ``asegstats2table`` command.
+    """
+    os.environ["SUBJECTS_DIR"] = subjects_dir
+    cmd = ['asegstats2table']
+
+    return cmd
+
+@wutils.cmdwrapper
+def aparcstats2table(subjects_dir, mask, tablefile, subjects, hemi, atlas,
+                     **kwargs):
+    """Wrapper for the ``aparcstats2table`` command.
+    """
+    os.environ["SUBJECTS_DIR"] = subjects_dir
+    cmd = ['aparcstats2table']
 
-    cmd += wutils.applyArgStyle('-', valmap=valmap, **kwargs)
- 
     return cmd

bip/main.py

@@ -19,6 +19,7 @@ import bip
 from bip.utils import setup_logging
 from bip.pipelines.struct_T1       import struct_T1
 from bip.pipelines.struct_T2_FLAIR import struct_T2_FLAIR
+from bip.pipelines.struct_FS       import struct_FS
 from bip.pipelines.struct_swMRI    import struct_swMRI
 from bip.pipelines.struct_asl      import struct_asl
 from bip.pipelines.dMRI_fieldmap   import dMRI_fieldmap
@@ -335,6 +336,7 @@ def main():
 
     pipe, targets = struct_T1.add_to_pipeline(ctx, pipe, tree, targets)
     pipe, targets = struct_T2_FLAIR.add_to_pipeline(ctx, pipe, tree, targets)
+    pipe, targets = struct_FS.add_to_pipeline(ctx, pipe, tree, targets)
     pipe, targets = struct_swMRI.add_to_pipeline(ctx, pipe, tree, targets)
     pipe, targets = struct_asl.add_to_pipeline(ctx, pipe, tree, targets)
     pipe, targets = dMRI_fieldmap.add_to_pipeline(ctx, pipe, tree, targets)

bip/pipelines/struct_FS/FS_get_IDPs.py

+#!/usr/bin/env python
+#
+# asl_get_IDPs.py - Sub-pipeline generating the IDPs of the ASL pipeline.
+#
+# 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
+# pylint: disable=W0613,R1718
+#
+
+import logging
+from pipe_tree import In, Out, Ref
+from bip.utils import redirect_logging
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        region_analysis_gm_csv:  In,
+        logs_dir:                Ref,
+        ASL_region_analysis_dir: Ref,
+        ASL_IDPs:                Out):
+
+    with redirect_logging('asl_get_IDPs', outdir=logs_dir):
+
+        result = ctx.subject
+
+        asl_format_file = ctx.get_data("asl/asl_format.txt")
+
+        with open(asl_format_file, 'rt', encoding="utf-8") as f:
+            asl_format = [x.strip().split() for x in f.readlines()]
+
+        regions = set([x[0] for x in asl_format])
+
+        for region in regions:
+            file_name_1 = ASL_region_analysis_dir + "/" + region + ".csv"
+            file_name_2 = ASL_region_analysis_dir + "/" + region + ".txt"
+
+            with open(file_name_1, 'rt', encoding="utf-8") as f:
+                file_data = f.read()
+
+            file_data = file_data.replace(", ", "_-_")
+            file_data = file_data.replace('"', "")
+            file_data = file_data.replace(" ", "_")
+
+            with open(file_name_2, 'wt', encoding="utf-8") as f:
+                f.write(file_data)
+
+
+        for elem in asl_format:
+
+            fil    = elem[0]
+            row    = elem[1]
+            column = elem[2]
+
+            file_name = ASL_region_analysis_dir + "/" + fil + ".txt"
+
+            with open(file_name, 'rt', encoding="utf-8") as f:
+                file_data = [x.strip().split(",") for x in f.readlines()]
+
+            for line in file_data:
+                if line[0] == row:
+                    val = line[int(column)]
+            if val == "":
+                val = "NaN"
+
+            result += " " + val
+
+        with open(ASL_IDPs, 'wt', encoding="utf-8") as f:
+            f.write(result + "\n")
+           
\ No newline at end of file

bip/pipelines/struct_FS/FS_proc.py

+#!/usr/bin/env python
+#
+# FS_proc.py - Sub-pipeline with the FreeSurfer main processing.
+#
+# 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
+# pylint: disable=W0613
+#
+
+import os
+import logging
+from pipe_tree import In, Out, Ref
+from bip import ext_wrappers
+from bip.utils import redirect_logging
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1_unbiased:                In,
+        T2_FLAIR_unbiased:          In(optional=True),
+        logs_dir:                   Ref,
+        rh_entorhinal_exvivo_label: Out):
+
+    with redirect_logging('FS_proc', outdir=logs_dir):
+        opt_T2_FLAIR = None
+
+        if os.path.exists(T2_FLAIR_unbiased):
+            opt_T2_FLAIR = T2_FLAIR_unbiased
+
+        ext_wrappers.recon_all(subjects_dir=os.getcwd(), directive="all",
+                               subjid="FreeSurfer", infile=T1_unbiased,
+                               FLAIR=opt_T2_FLAIR)
+        print("PATATA")

bip/pipelines/struct_FS/FS_segm.py

+#!/usr/bin/env python
+#
+# asl_proc.py - Sub-pipeline with the ASL main processing.
+#
+# 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
+# pylint: disable=W0613
+#
+
+import os
+import glob
+import logging
+from shutil import copyfile
+from fsl import wrappers
+from pipe_tree import In, Out, Ref
+from gradunwarp.core.gradient_unwarp_apply import gradient_unwarp_apply
+from bip.utils import redirect_logging
+from bip.commands.bb_read_json_field import bb_read_json_field
+
+log = logging.getLogger(__name__)
+
+def run(ctx,
+        T1:                            In,
+        T1_brain:                      In,
+        T1_to_MNI_warp:                In,
+        fieldmap_fout_to_T1_brain_rad: In,
+        ASL_M0:                        In,
+        ASL_M0_json:                   In,
+        logs_dir:                      Ref,
+        ASL_GDC:                       Ref,
+        OXASL_ra_dir:                  Ref,
+        T1_fast_prefix:                Ref,
+        ASL_PLD_prefix:                Ref,
+        ASL_raw_dir:                   Ref,
+        ASL_M0_ud:                     Out,
+        ASL_M0_ud_warp:                Out,
+        ASL_DATA_wrongorder:           Out,
+        ASL_DATA:                      Out,
+        ASL_DATA_diff:                 Out,
+        ASL_DATA_diff_mean:            Out,
+        ASL_control:                   Out,
+        ASL_label:                     Out,
+        CALIB:                         Out,
+        CALIB_json:                    Out,
+        ASL_qc_modelfit:               Out,
+        ASL_mni_right_amygdala:        Out,
+        ASL_masks_right_amygdala:      Out,
+        ASL_native_right_amygdala:     Out,
+        ASL_struct_90_wm:              Out,
+        ASL_calib_M0:                  Out,
+        ASL_struct_ACBV_calib:         Out,
+        ASL_std_ACBV_calib:            Out,
+        region_analysis_gm_csv:        Out):
+
+    with redirect_logging('asl_proc', outdir=logs_dir):
+
+        BBASL_ROI_DIR = ctx.get_data("asl/ukb_rois/")
+
+        ###############################################
+        # Initial check for Phase Encoding Directions #
+        ###############################################
+        list_PEDs = []
+
+        for fil in glob.glob(ASL_raw_dir + "/*.json"):
+            list_PEDs.append(bb_read_json_field(fileName=fil,
+                                fieldName="PhaseEncodingDirection"))
+
+        if len(list_PEDs) != 12:
+            error_msg =  "ERROR: Wrong phase-encoding-direction in some of the "
+            error_msg += "files (or wrong number of files). "
+            error_msg += "ASL will not be processed."
+            print(error_msg)
+            # TODO: Make unusable this modality
+            # bb_make_unusable ASL $src "2 Wrong_acquisition"
+
+        ##################
+        # Pre-processing #
+        ##################
+        list_fils_control = glob.glob(ASL_PLD_prefix + "*_control.nii.gz")
+        list_fils_control.sort()
+        list_fils_label = glob.glob(ASL_PLD_prefix + "*_label.nii.gz")
+        list_fils_label.sort()
+
+        wrappers.fslmerge("t", ASL_control, *list_fils_control)
+        wrappers.fslmerge("t", ASL_label, *list_fils_label)
+        wrappers.fslmerge("t", ASL_DATA_wrongorder, ASL_label, ASL_control)
+
+        os.symlink("../raw/ASL_M0.nii.gz", CALIB)
+        os.symlink("../raw/ASL_M0.json",   CALIB_json)
+
+        TE = bb_read_json_field(fileName=CALIB_json, fieldName="EchoTime",
+                                rounding = 3, multFactor=1000)
+
+        #Gradient distortion correction applied to the M0
+        if ctx.gdc != '':
+            #Calculate and apply the Gradient Distortion Unwarp
+            # TODO: Review the "half=True" in next version
+            gradient_unwarp_apply(WD=ASL_GDC, infile=CALIB, outfile=ASL_M0_ud,
+                                  owarp=ASL_M0_ud_warp,gradcoeff=ctx.gdc,
+                                  vendor='siemens', nojac=True, half=False)
+        else:
+            copyfile(src=CALIB, dst=ASL_M0_ud)
+
+
+        ##############
+        # Processing #
+        ##############
+        # Use asl_file to get into alternating tag-control
+        wrappers.asl_file(data=ASL_DATA_wrongorder, out=ASL_DATA, ntis=5,
+                          iaf="tcb")
+
+        # PWI generation - this is more of a sanity check than anything else
+        # now use ASL file to get difference images
+        wrappers.asl_file(data=ASL_DATA, out=ASL_DATA_diff,
+                          mean=ASL_DATA_diff_mean, ntis=5, iaf="tc", diff=True)
+
+
+        #######################
+        # PERFUSION ANALSYSIS #
+        #######################
+
+        raa_list = [BBASL_ROI_DIR + "MNI_seg_max_prob_masked_RandL.nii.gz",
+                    BBASL_ROI_DIR + "HO_L_Cerebral_WM_thr80.nii.gz",
+                    BBASL_ROI_DIR + "HO_R_Cerebral_WM_thr80.nii.gz",
+                    BBASL_ROI_DIR + "VascularTerritories_ero.nii.gz"]
+        raa_labels_list = [BBASL_ROI_DIR + "MNI_seg_max_prob_masked_RandL.txt",
+                           BBASL_ROI_DIR + "HO_L_Cerebral_WM_thr80.txt",
+                           BBASL_ROI_DIR + "HO_R_Cerebral_WM_thr80.txt",
+                           BBASL_ROI_DIR + "VascularTerritories_ero.txt"]
+
+        wrappers.oxford_asl(data=ASL_DATA,
+                            out=OXASL_ra_dir,
+                            c=CALIB,
+                            s=T1,
+                            regfrom_method="pwi",
+                            reg_init_bbr=True,
+                            sbrain=T1_brain,
+                            fastsrc=T1_fast_prefix,
+                            warp=T1_to_MNI_warp,
+                            fmap=fieldmap_fout_to_T1_brain_rad,
+                            fmapmag=T1,
+                            fmapmagbrain=T1_brain,
+                            gdcwarp=ASL_M0_ud_warp,
+                            iaf="tc",
+                            ibf="rpt",
+                            mc="on",
+                            tis=[2.2, 2.6, 3.0, 3.4, 3.8],
+                            casl=True,
+                            bolus=1.8,
+                            fixbolus=True,
+                            cgain=10,
+                            spatial=True,
+                            tr=5,
+                            te=TE,
+                            echospacing=0.0005,
+                            pedir="x",
+                            gm_thresh=0.7,
+                            cmethod="voxel",
+                            nofmapreg=True,
+                            region_analysis=True,
+                            region_analysis_atlas= raa_list,
+                            region_analysis_atlas_labels= raa_labels_list,
+                            region_analysis_save_rois=True,
+                            qc_output=True)

bip/pipelines/struct_FS/struct_FS.py

+#!/usr/bin/env python
+#
+# struct_FS.py - Pipeline with the FreeSurfer processing.
+#
+# 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 logging
+from bip.utils import redirect_logging
+from bip.pipelines.struct_FS import FS_proc#, FS_segm, FS_get_IDPs
+
+log = logging.getLogger(__name__)
+
+def add_to_pipeline(ctx, pipe, tree, targets):
+
+    logs_dir=tree.get('logs_dir')
+
+    with redirect_logging('pipe_struct_FS', outdir=logs_dir):
+        pipe(FS_proc.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        targets.append('rh_entorhinal_exvivo_label')
+        #pipe(FS_segm.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        #targets.append('ASL_IDPs')
+        #pipe(FS_get_IDPs.run, submit=dict(jobtime=200), kwargs={'ctx' : ctx})
+        #targets.append('FS_IDPs')
+
+    return pipe, targets

bip/pipelines/struct_asl/struct_asl.py

@@ -11,8 +11,7 @@
 
 import logging
 from bip.utils import redirect_logging
-from bip.pipelines.struct_asl import asl_proc
-from bip.pipelines.struct_asl import asl_get_IDPs
+from bip.pipelines.struct_asl import asl_proc, asl_get_IDPs
 
 log = logging.getLogger(__name__)
 

init_vars

 #!/bin/bash
 
+# For the pipeline. Requires the existence of a SWdir with:
+#  - The pipeline
+#  - fsl
+#  - freesurfer 
 export BB_BIN_DIR="$SWdir/bip"
 export FSLDIR="$SWdir/fsl"
 
+# For FSL
 source $FSLDIR/etc/fslconf/fsl.sh
-source $FSLDIR/bin/activate ukb
 
 export FSLOUTPUTTYPE="NIFTI_GZ"
 export FSLSUB_CONF="$SWdir/bip/fsl_sub.yml"
 
+# For FreeSurfer
+unset FSFAST_HOME
+unset MNI_DIR
+unset FSL_DIR
+export FREESURFER_HOME="$SWdir/freesurfer_7.3.2"
+export MNI_DIR="${FREESURFER_HOME}/mni"
+source $FREESURFER_HOME/SetUpFreeSurfer.sh
+
+# Activate the environment
+source $FSLDIR/bin/activate ukb
+
+# Unsetting some problematic variables
 unset PYTHONPATH
 unset LC_ALL
 unset LC_CTYPE
 
+# Creating a useful alias
 alias runbip="$BB_BIN_DIR/bip/main.py"
 PS1="(ukb) \! "`whoami`@`hostname|cut -d. -f1`@$(basename $(tty))" | \t | \w $ " ; 
 export PS1
\ No newline at end of file