Merge pull request #53 from akeshavan/incorporate_eddy_quad

Incorporate eddy quad

Author: akeshavan

.gitignore

@@ -103,3 +103,6 @@ ENV/
 
 # Mac OS nonsense:
 .DS_Store
+
+#kubernetes stuff
+kubernetes/jobs/

dmriprep/data.py

@@ -425,6 +425,21 @@ def postprocess(self, subject):
                             stdout=fnull,
                             stderr=subprocess.STDOUT)
 
+            # if the site is CBIC, then the freesurfer directory has an additional level.
+            # move that level up by 1 (e.g. removing the T1w_VNavNorm folder
+            if subject.site == 'Site-CBIC':
+                newpath = freesurfer_path.replace('T1w_VNavNorm/', '')
+                move_cmd = 'mv {oldpath} {newpath}'.format(oldpath=freesurfer_path, newpath=newpath)
+                fnull1 = open(os.devnull, 'w')
+                subprocess.call(move_cmd.split(),
+                                stdout=fnull1,
+                                stderr=subprocess.STDOUT)
+
+            # now check that the AP/PA files are named correctly
+            # eg it should look like "sub-{id}_dir-{dir}_acq-dwi_epi.nii.gz
+            # but sometimes it looks like sub-{id}_acq-dwi_run-01_epi.nii.gz
+            # which is silly. the direction should be in the filename.
+
 
 class Subject:
     """A single dMRI study subject"""

dmriprep/qc.py

@@ -3,7 +3,7 @@
 import os.path as op
 from dipy.segment.mask import median_otsu
 from io import BytesIO
-from nipype.utils.filemanip import save_json
+from nipype.utils.filemanip import save_json, load_json
 import base64
 import matplotlib
 matplotlib.use('agg')
@@ -212,6 +212,7 @@ def createB0_ColorFA_Mask_Sprites(b0_file, colorFA_file, mask_file):
 def create_report_json(dwi_corrected_file, eddy_rms, eddy_report,
                        color_fa_file, anat_mask_file,
                        outlier_indices,
+                       eddy_qc_file,
                        outpath=op.abspath('./report.json')):
 
     report = {}
@@ -229,6 +230,7 @@ def create_report_json(dwi_corrected_file, eddy_rms, eddy_report,
         report['eddy_report'] = f.readlines()
 
     report['eddy_params'] = np.genfromtxt(eddy_rms).tolist()
-
+    eddy_qc = load_json(eddy_qc_file)
+    report['eddy_quad'] = eddy_qc
     save_json(outpath, report)
     return outpath

dmriprep/run.py

@@ -301,7 +301,8 @@ def run_dmriprep_pe(subject_id, dwi_file, dwi_file_AP, dwi_file_PA,
 
     # write the graph (this is saved to the working dir)
     wf.write_graph()
-
+    wf.config['execution']['remove_unnecessary_outputs'] = False
+    wf.config['execution']['keep_inputs'] = True
     wf.run()
 
 
@@ -363,6 +364,31 @@ def get_dmriprep_pe_workflow():
     import multiprocessing
     eddy.inputs.num_threads = multiprocessing.cpu_count()
 
+    eddy_quad = pe.Node(fsl.EddyQuad(verbose=True), name="eddy_quad")
+    get_path = lambda x: x.split('.nii.gz')[0]
+    wf.connect(prep, ('fsl_eddy.out_corrected', get_path), eddy_quad, "base_name")
+    wf.connect(inputspec, 'bval_file', eddy_quad, 'bval_file')
+    wf.connect(prep, 'Rotate_Bvec.out_file', eddy_quad, 'bvec_file')
+    wf.connect(prep, 'peb_correction.topup.out_field', eddy_quad, 'field')
+    wf.connect(prep, 'gen_index.out_file', eddy_quad, 'idx_file')
+    wf.connect(prep, 'peb_correction.topup.out_enc_file', eddy_quad, 'param_file')
+
+    # need a mask file for eddy_quad. lets get it from the B0.
+    def get_b0_mask_fn(b0_file):
+        import nibabel as nib
+        from nipype.utils.filemanip import fname_presuffix
+        from dipy.segment.mask import median_otsu
+        import os
+
+        mask_file = fname_presuffix(b0_file, suffix="_mask", newpath=os.path.abspath('.'))
+        img = nib.load(b0_file)
+        data, aff = img.get_data(), img.affine
+        _, mask = median_otsu(data, 2, 1)
+        nib.Nifti1Image(mask.astype(float), aff).to_filename(mask_file)
+        return mask_file
+
+
+
     def id_outliers_fn(outlier_report, threshold, dwi_file):
         """Get list of scans that exceed threshold for number of outliers
 
@@ -437,13 +463,19 @@ def num_outliers(scan, outliers):
     wf.connect(inputspec, 'dwi_file_pa', list_merge, 'in2')
 
     merge = pe.Node(fsl.Merge(dimension='t'), name="mergeAPPA")
-    # merge.inputs.in_files = [dwi_file_ap, dwi_file_pa]
     wf.connect(merge, 'merged_file', prep, 'inputnode.alt_file')
     wf.connect(list_merge, 'out', merge, 'in_files')
 
     fslroi = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name="fslroi")
     wf.connect(prep, "outputnode.out_file", fslroi, "in_file")
 
+    b0mask_node = pe.Node(niu.Function(input_names=['b0_file'],
+                                       output_names=['mask_file'],
+                                       function=get_b0_mask_fn),
+                          name="getB0Mask")
+    wf.connect(fslroi, 'roi_file', b0mask_node, 'b0_file')
+    wf.connect(b0mask_node, 'mask_file', eddy_quad, 'mask_file')
+
     bbreg = pe.Node(fs.BBRegister(contrast_type="t2", init="coreg",
                                   out_fsl_file=True,
                                   # subjects_dir=subjects_dir,
@@ -487,7 +519,7 @@ def drop_outliers_fn(in_file, in_bval, in_bvec, drop_scans):
         from nipype.utils.filemanip import fname_presuffix
 
         img = nib.load(op.abspath(in_file))
-        img_data = img.get_fdata()
+        img_data = img.get_data()
         img_data_thinned = np.delete(img_data,
                                      drop_scans,
                                      axis=3)
@@ -679,7 +711,7 @@ def binarize_aparc(aparc_aseg):
     wf.connect(prep, "fsl_eddy.out_residuals",
                datasink, "dmriprep.qc.@eddyresid")
 
-    # the file that told us which volumes to trop
+    # the file that told us which volumes to drop
     wf.connect(id_outliers_node, "outpath", datasink, "dmriprep.qc.@droppedscans")
 
     # the tensors of the dropped volumes dwi
@@ -700,6 +732,15 @@ def binarize_aparc(aparc_aseg):
     wf.connect(get_tensor_eddy, "color_fa_file", datasink, "dmriprep.dti_eddy.@colorfa")
     wf.connect(scale_tensor_eddy, "out_file", datasink, "dmriprep.dti_eddy.@scaled_tensor")
 
+    # all the eddy_quad stuff
+    wf.connect(eddy_quad, 'out_qc_json', datasink, "dmriprep.qc.@eddyquad_json")
+    wf.connect(eddy_quad, 'out_qc_pdf', datasink, "dmriprep.qc.@eddyquad_pdf")
+    wf.connect(eddy_quad, 'out_avg_b_png', datasink, "dmriprep.qc.@eddyquad_bpng")
+    wf.connect(eddy_quad, 'out_avg_b0_png', datasink, "dmriprep.qc.@eddyquad_b0png")
+    wf.connect(eddy_quad, 'out_cnr_png', datasink, "dmriprep.qc.@eddyquad_cnr")
+    wf.connect(eddy_quad, 'out_vdm_png', datasink, "dmriprep.qc.@eddyquad_vdm")
+    wf.connect(eddy_quad, 'out_residuals', datasink, 'dmriprep.qc.@eddyquad_resid')
+
     # anatomical registration stuff
     wf.connect(bbreg, "min_cost_file", datasink, "dmriprep.reg.@mincost")
     wf.connect(bbreg, "out_fsl_file", datasink, "dmriprep.reg.@fslfile")
@@ -711,23 +752,26 @@ def binarize_aparc(aparc_aseg):
     wf.connect(reslice_orig_to_dwi, 'out_file', datasink, 'dmriprep.anat.@T1w')
 
     def report_fn(dwi_corrected_file, eddy_rms, eddy_report,
-                  color_fa_file, anat_mask_file, outlier_indices):
+                  color_fa_file, anat_mask_file, outlier_indices,
+                  eddy_qc_file):
         from dmriprep.qc import create_report_json
 
         report = create_report_json(dwi_corrected_file, eddy_rms, eddy_report,
-                                    color_fa_file, anat_mask_file, outlier_indices)
+                                    color_fa_file, anat_mask_file, outlier_indices,
+                                    eddy_qc_file)
         return report
 
     report_node = pe.Node(niu.Function(
         input_names=['dwi_corrected_file', 'eddy_rms',
                      'eddy_report', 'color_fa_file',
-                     'anat_mask_file', 'outlier_indices'],
+                     'anat_mask_file', 'outlier_indices', 'eddy_qc_file'],
         output_names=['report'],
         function=report_fn
     ), name="reportJSON")
 
     # for the report, lets show the eddy corrected (full volume) image
     wf.connect(voltransform, "transformed_file", report_node, 'dwi_corrected_file')
+    wf.connect(eddy_quad, 'out_qc_json', report_node, 'eddy_qc_file')
 
     # add the rms movement output from eddy
     wf.connect(prep, "fsl_eddy.out_movement_rms", report_node, 'eddy_rms')
@@ -762,7 +806,6 @@ def name_files_nicely(dwi_file, subject_id):
             ("art.eddy_corrected_outliers", dwi_fname.replace("dwi", "outliers")),
             ("color_fa", "colorfa"),
             ("orig_out", dwi_fname.replace("_dwi", "_T1w")),
-            # ("eddy_corrected_", dwi_fname.replace("dwi", "")),
             ("stats.eddy_corrected", dwi_fname.replace("dwi", "artStats")),
             ("eddy_corrected.eddy_parameters", dwi_fname + ".eddy_parameters"),
             ("qc/eddy_corrected", "qc/" + dwi_fname),

docker/Dockerfile

@@ -39,8 +39,8 @@ RUN export ND_ENTRYPOINT="/neurodocker/startup.sh" \
 
 ENTRYPOINT ["/neurodocker/startup.sh"]
 
-ENV FSLDIR="/opt/fsl-6.0.0" \
-    PATH="/opt/fsl-6.0.0/bin:$PATH"
+ENV FSLDIR="/opt/fsl-5.0.11" \
+    PATH="/opt/fsl-5.0.11/bin:$PATH"
 RUN apt-get update -qq \
     && apt-get install -y -q --no-install-recommends \
            bc \
@@ -63,15 +63,15 @@ RUN apt-get update -qq \
     && apt-get clean \
     && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* \
     && echo "Downloading FSL ..." \
-    && mkdir -p /opt/fsl-6.0.0 \
-    && curl -fsSL --retry 5 https://fsl.fmrib.ox.ac.uk/fsldownloads/fsl-6.0.0-centos6_64.tar.gz \
-    | tar -xz -C /opt/fsl-6.0.0 --strip-components 1 \
+    && mkdir -p /opt/fsl-5.0.11 \
+    && curl -fsSL --retry 5 https://fsl.fmrib.ox.ac.uk/fsldownloads/fsl-5.0.11-centos6_64.tar.gz \
+    | tar -xz -C /opt/fsl-5.0.11 --strip-components 1 \
     && sed -i '$iecho Some packages in this Docker container are non-free' $ND_ENTRYPOINT \
     && sed -i '$iecho If you are considering commercial use of this container, please consult the relevant license:' $ND_ENTRYPOINT \
     && sed -i '$iecho https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Licence' $ND_ENTRYPOINT \
     && sed -i '$isource $FSLDIR/etc/fslconf/fsl.sh' $ND_ENTRYPOINT \
     && echo "Installing FSL conda environment ..." \
-    && bash /opt/fsl-6.0.0/etc/fslconf/fslpython_install.sh -f /opt/fsl-6.0.0
+    && bash /opt/fsl-5.0.11/etc/fslconf/fslpython_install.sh -f /opt/fsl-5.0.11
 
 ENV CONDA_DIR="/opt/miniconda-latest" \
     PATH="/opt/miniconda-latest/bin:$PATH"
@@ -134,8 +134,8 @@ COPY ./license.txt /opt/freesurfer-6.0.0/license.txt
 
 #&& sync && conda clean -tipsy && sync
 
-ADD environment.yml environment.yml
 RUN apt-get update && apt-get install -y git gcc libopenblas-base
+ADD environment.yml environment.yml
 RUN conda env create -f environment.yml
 
 ENV LD_LIBRARY_PATH=/usr/lib/openblas-base/
@@ -153,7 +153,7 @@ RUN echo '{ \
     \n    [ \
     \n      "fsl", \
     \n      { \
-    \n        "version": "6.0.0" \
+    \n        "version": "5.0.11" \
     \n      } \
     \n    ], \
     \n    [ \

docker/environment.yml

@@ -11,6 +11,7 @@ dependencies:
    - pandas
    - tqdm
    - pip:
-     - "--editable=git+https://github.com/nipy/nipype@35b58cdb59a900430f4b9639e220c8408341805d#egg=nipype"
+     - "--editable=git+https://github.com/akeshavan/nipype@7b5e7ab92183eb561980ec713d86eb69f28733e1#egg=nipype"
+     - "--editable=git+https://git.fmrib.ox.ac.uk/matteob/eddy_qc_release@57bb11da6a634c4195593fbc439ba9f8998157b0#egg=eddy_qc"
      - bids
      - duecredit

kubernetes/create_kube_job.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+import sys
+from boto3 import Session
+
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        subject = sys.argv[1]
+
+        session = Session()
+        credentials = session.get_credentials()
+        # Credentials are refreshable, so accessing your access key / secret key
+        # separately can lead to a race condition. Use this to get an actual matched
+        # set.
+        current_credentials = credentials.get_frozen_credentials()
+
+        access_key = current_credentials.access_key
+        secret_key = current_credentials.secret_key
+
+        with open("run_dmriprep.yml.tmpl", 'r') as template:
+            with open("jobs/job_{}.yml".format(subject), 'w') as f:
+                all_text = "\n".join(template.readlines())
+                all_text = all_text.replace("{{subject_lower}}", subject.lower())
+                all_text = all_text.replace("{{subject}}", subject)
+                all_text = all_text.replace("{{access_key}}", access_key)
+                all_text = all_text.replace("{{secret_key}}", secret_key)
+                f.write(all_text)
+

kubernetes/docker/build_tag_push.sh

@@ -1,4 +1,4 @@
 #!/bin/bash
-docker build -t dmriprep:kube1 -f dockerfile-dmriprep-kube1 .
-docker tag dmriprep:kube1 gcr.io/dmriprep/dmriprep:kube1
-docker push gcr.io/dmriprep/dmriprep:kube1
+docker build -t dmriprep:kube3 -f dockerfile-dmriprep-kube .
+docker tag dmriprep:kube3 gcr.io/dmriprep/dmriprep:kube3
+docker push gcr.io/dmriprep/dmriprep:kube3

kubernetes/run_dmriprep.yml renamed to kubernetes/run_dmriprep.yml.tmpl

@@ -1,13 +1,13 @@
 apiVersion: batch/v1
 kind: Job
 metadata:
-  name: dmriprep0
+  name: dmriprep-{{subject_lower}}
 spec:
   template:
     spec:
       containers:
       - name: dmriprep-subject
-        image: gcr.io/dmriprep/dmriprep:kube1
-        command: ["/neurodocker/startup.sh", "./dmriprep_all.sh",  "sub-NDARED632KNG", "{{access_key}}", "{{secret_key}}"]
+        image: gcr.io/dmriprep/dmriprep:kube3
+        command: ["/neurodocker/startup.sh", "./dmriprep_all.sh",  "{{subject}}", "{{access_key}}", "{{secret_key}}"]
       restartPolicy: Never
   backoffLimit: 1

kubernetes/setup_gcp_kubernetes_dmriprep.sh

@@ -15,4 +15,4 @@ CLUSTERNAME=dmriprep
 gcloud config set compute/zone $ZONE
 
 # start the cluster!
-gcloud beta container clusters create $CLUSTERNAME --machine-type n1-highmem-4 --enable-autoscaling --max-nodes=$MAX_NODES --num-nodes 2 --cluster-version latest --node-labels dmriprep/node-purpose=core
+gcloud beta container clusters create $CLUSTERNAME --machine-type n1-highmem-4 --enable-autoscaling --max-nodes=$MAX_NODES --num-nodes 1 --cluster-version latest --node-labels dmriprep/node-purpose=core