Merge pull request #38 from fliem/euler

Euler number
bids-apps · Dec 1, 2017 · cf36097 · cf36097
2 parents 81fde04 + 860fc8b
commit cf36097
Show file tree

Hide file tree

Showing 4 changed files with 96 additions and 22 deletions.
diff --git a/Dockerfile b/Dockerfile
@@ -19,7 +19,7 @@ RUN wget -qO- https://surfer.nmr.mgh.harvard.edu/pub/dist/freesurfer/6.0.1/frees
 
 RUN apt-get install -y python3
 RUN apt-get install -y python3-pip
-RUN pip3 install nibabel
+RUN pip3 install nibabel pandas
 RUN apt-get install -y python2.7
 RUN apt-get install -y python-pip
 

diff --git a/README.md b/README.md
@@ -25,11 +25,15 @@ This App has the following command line arguments:
         $ docker run -ti --rm bids/freesurfer --help
         usage: run.py [-h]
                       [--participant_label PARTICIPANT_LABEL [PARTICIPANT_LABEL ...]]
+                      [--session_label SESSION_LABEL [SESSION_LABEL ...]]
                       [--n_cpus N_CPUS]
-                      [--stages {autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon2-pial,autorecon3,autorecon-all,all}
-                                [{autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon2-pial,autorecon3,autorecon-all,all} ...]]
-                      [--template_name TEMPLATE_NAME] --license_key LICENSE_KEY
+                      [--stages {autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon-pial,autorecon3,autorecon-all,all}
+                                [{autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon-pial,autorecon3,autorecon-all,all} ...]]
+                      [--steps {cross-sectional,template,longitudinal}
+                               [{cross-sectional,template,longitudinal} ...]]
+                      [--template_name TEMPLATE_NAME] --license_file LICENSE_FILE
                       [--acquisition_label ACQUISITION_LABEL]
+                      [--refine_pial_acquisition_label REFINE_PIAL_ACQUISITION_LABEL]
                       [--multiple_sessions {longitudinal,multiday}]
                       [--refine_pial {T2,FLAIR,None,T1only}]
                       [--hires_mode {auto,enable,disable}]
@@ -39,6 +43,7 @@ This App has the following command line arguments:
                       [-v] [--bids_validator_config BIDS_VALIDATOR_CONFIG]
                       [--skip_bids_validator]
                       bids_dir output_dir {participant,group1,group2}
+
         FreeSurfer recon-all + custom template generation.
 
         positional arguments:
@@ -51,10 +56,10 @@ This App has the following command line arguments:
           {participant,group1,group2}
                                 Level of the analysis that will be performed. Multiple
                                 participant level analyses can be run independently
-                                (in parallel) using the same output_dir. "goup1"
-                                creates study specific group template. "group2 exports
-                                group stats tables for cortical parcellation and
-                                subcortical segmentation.
+                                (in parallel) using the same output_dir. "group1"
+                                creates study specific group template. "group2"
+                                exports group stats tables for cortical parcellation,
+                                subcortical segmentation a table with euler numbers.
 
         optional arguments:
           -h, --help            show this help message and exit
@@ -65,21 +70,34 @@ This App has the following command line arguments:
                                 parameter is not provided all subjects should be
                                 analyzed. Multiple participants can be specified with
                                 a space separated list.
+          --session_label SESSION_LABEL [SESSION_LABEL ...]
+                                The label of the session that should be analyzed. The
+                                label corresponds to ses-<session_label> from the BIDS
+                                spec (so it does not include "ses-"). If this
+                                parameter is not provided all sessions should be
+                                analyzed. Multiple sessions can be specified with a
+                                space separated list.
           --n_cpus N_CPUS       Number of CPUs/cores available to use.
-          --stages {autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon2-pial,autorecon3,autorecon-all,all}
-                                [{autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon2-pial,autorecon3,autorecon-all,all} ...]
+          --stages {autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon-pial,autorecon3,autorecon-all,all}
+                                [{autorecon1,autorecon2,autorecon2-cp,autorecon2-wm,autorecon-pial,autorecon3,autorecon-all,all} ...]
                                 Autorecon stages to run.
+          --steps {cross-sectional,template,longitudinal} [{cross-sectional,template,longitudinal} ...]
+                                Longitudinal pipeline steps to run.
           --template_name TEMPLATE_NAME
                                 Name for the custom group level template generated for
                                 this dataset
           --license_file LICENSE_FILE
                                 Path to FreeSurfer license key file. To obtain it you
-								need to register (for free) visit
+                                need to register (for free) at
                                 https://surfer.nmr.mgh.harvard.edu/registration.html
           --acquisition_label ACQUISITION_LABEL
                                 If the dataset contains multiple T1 weighted images
                                 from different acquisitions which one should be used?
                                 Corresponds to "acq-<acquisition_label>"
+          --refine_pial_acquisition_label REFINE_PIAL_ACQUISITION_LABEL
+                                If the dataset contains multiple T2 or FLAIR weighted
+                                images from different acquisitions which one should be
+                                used? Corresponds to "acq-<acquisition_label>"
           --multiple_sessions {longitudinal,multiday}
                                 For datasets with multiday sessions where you do not
                                 want to use the longitudinal pipeline, i.e., sessions
@@ -100,15 +118,15 @@ This App has the following command line arguments:
                                 stats from.
           --measurements {area,volume,thickness,thicknessstd,meancurv,gauscurv,foldind,curvind}
                                 [{area,volume,thickness,thicknessstd,meancurv,gauscurv,foldind,curvind} ...]
-                                Group2 option: cortical measurements to extract stats for.
+                                Group2 option: cortical measurements to extract stats
+                                for.
           -v, --version         show program's version number and exit
           --bids_validator_config BIDS_VALIDATOR_CONFIG
                                 JSON file specifying configuration of bids-validator:
                                 See https://github.com/INCF/bids-validator for more
                                 info
           --skip_bids_validator
-                               skips bids validation
-
+                                skips bids validation
 
 #### Participant level
 To run it in participant level mode (for one participant):
@@ -125,6 +143,7 @@ To run it in participant level mode (for one participant):
 After doing this for all subjects (potentially in parallel) the
 group level analyses can be run.
 
+##### Template creation
 To create a study specific template run:
 
 		docker run -ti --rm \
@@ -134,13 +153,24 @@ To create a study specific template run:
 		/bids_dataset /outputs group1 \
 		--license_file "license.txt"
 
+##### Stats and quality tables export
 To export tables with aggregated measurements within regions of
-cortical parcellation and subcortical segementation run:
+cortical parcellation and subcortical segementation, and a table with
+ euler numbers (a quality metric, see
+ [Rosen et. al, 2017](https://www.biorxiv.org/content/early/2017/10/01/125161))
+ run:
 
 		docker run -ti --rm \
 		-v /Users/filo/data/ds005:/bids_dataset:ro \
 		-v /Users/filo/outputs:/outputs \
 		bids/freesurfer \
 		/bids_dataset /outputs group2 \
 		--license_file "license.txt"
-Also see *--parcellations* and *--measurements* arguments.
+Also see the *--parcellations* and *--measurements* arguments.
+
+This step writes ouput into `<output_dir>/00_group2_stats_tables/`. E.g.:
+
+* `lh.aparc.thickness.tsv` contains cortical thickness values for the
+left hemisphere extracted via the aparac parcellation.
+* `aseg.tsv` contains subcortical information from the aseg segmentation.
+* `euler.tsv` contains the euler numbers
diff --git a/circle.yml b/circle.yml
@@ -28,14 +28,15 @@ test:
   override:
     # print version
     - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test1:/bids_dataset bids/${CIRCLE_PROJECT_REPONAME,,} --version
+    - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test1:/bids_dataset bids/${CIRCLE_PROJECT_REPONAME,,} -h
     - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test1:/bids_dataset -v ${HOME}/outputs1:/outputs bids/${CIRCLE_PROJECT_REPONAME,,} /bids_dataset /outputs participant --participant_label 01 --license_file=/license.txt --stages autorecon1 && cat ${HOME}/outputs1/sub-01/scripts/recon-all.done :
         timeout: 21600
     - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test2:/bids_dataset -v ${HOME}/outputs2:/outputs bids/${CIRCLE_PROJECT_REPONAME,,} /bids_dataset /outputs participant --participant_label 01 --steps cross-sectional --session_label test --license_file=/license.txt --stages autorecon1 && cat ${HOME}/outputs2/sub-01_ses-test/scripts/recon-all.done :
         timeout: 21600
     # group2 tests
-    - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test1:/bids_dataset -v ${HOME}/data/ds114_test1_freesurfer_precomp_v6.0.0:/outputs bids/${CIRCLE_PROJECT_REPONAME,,} /bids_dataset /outputs group2 --license_file=/license.txt  && mkdir -p ${HOME}/outputs1/ && sudo mv ${HOME}/data/ds114_test1_freesurfer_precomp_v6.0.0/00_group* ${HOME}/outputs1/ && cat ${HOME}/outputs1/00_group2_stats_tables/lh.aparc.thickness.tsv :
+    - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test1:/bids_dataset -v ${HOME}/data/ds114_test1_freesurfer_precomp_v6.0.0:/outputs bids/${CIRCLE_PROJECT_REPONAME,,} /bids_dataset /outputs group2 --license_file=/license.txt  && mkdir -p ${HOME}/outputs1/ && sudo mv ${HOME}/data/ds114_test1_freesurfer_precomp_v6.0.0/00_group* ${HOME}/outputs1/ && cat ${HOME}/outputs1/00_group2_stats_tables/lh.aparc.thickness.tsv && cat ${HOME}/outputs1/00_group2_stats_tables/euler.tsv:
         timeout: 21600
-    - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test2:/bids_dataset -v ${HOME}/data/ds114_test2_freesurfer_precomp_v6.0.0:/outputs bids/${CIRCLE_PROJECT_REPONAME,,} /bids_dataset /outputs group2 --license_file=/license.txt  && mkdir -p ${HOME}/outputs2/ && sudo mv ${HOME}/data/ds114_test2_freesurfer_precomp_v6.0.0/00_group* ${HOME}/outputs2/ && cat ${HOME}/outputs2/00_group2_stats_tables/lh.aparc.thickness.tsv :
+    - docker run -ti --rm --read-only -v $PWD/license.txt:/license.txt -v /tmp:/tmp -v /var/tmp:/var/tmp -v ${HOME}/data/ds114_test2:/bids_dataset -v ${HOME}/data/ds114_test2_freesurfer_precomp_v6.0.0:/outputs bids/${CIRCLE_PROJECT_REPONAME,,} /bids_dataset /outputs group2 --license_file=/license.txt  && mkdir -p ${HOME}/outputs2/ && sudo mv ${HOME}/data/ds114_test2_freesurfer_precomp_v6.0.0/00_group* ${HOME}/outputs2/ && cat ${HOME}/outputs2/00_group2_stats_tables/lh.aparc.thickness.tsv && cat ${HOME}/outputs2/00_group2_stats_tables/euler.tsv:
         timeout: 21600
 
 deployment:

diff --git a/run.py b/run.py
@@ -8,7 +8,8 @@
 from shutil import rmtree
 import subprocess
 from warnings import warn
-
+import pandas as pd
+import re
 
 def run(command, env={}, ignore_errors=False):
     merged_env = os.environ
@@ -39,7 +40,8 @@ def run(command, env={}, ignore_errors=False):
                     'Multiple participant level analyses can be run independently '
                     '(in parallel) using the same output_dir. '
                     '"group1" creates study specific group template. '
-                    '"group2 exports group stats tables for cortical parcellation and subcortical segmentation.',
+                    '"group2" exports group stats tables for cortical parcellation, subcortical segmentation '
+                                           'a table with euler numbers.',
                     choices=['participant', 'group1', 'group2'])
 parser.add_argument('--participant_label', help='The label of the participant that should be analyzed. The label '
                     'corresponds to sub-<participant_label> from the BIDS spec '
@@ -163,7 +165,6 @@ def run(command, env={}, ignore_errors=False):
     env = {'FS_LICENSE': args.license_file}
 else:
     raise Exception("Provided license file does not exist")
-    raise Exception("Provided license file does not exist")
 
 # running participant level
 if args.analysis_level == "participant":
@@ -491,4 +492,46 @@ def run(command, env={}, ignore_errors=False):
         print("\nTable export finished for %d subjects/sessions." % len(subjects))
 
     else:
-        print("\nNo subjects included in the analysis. Skipping group2 level.")
+        print("\nNo subjects included in the analysis. Skipping group2 level stats tables.")
+
+
+    # This extracts the euler numbers for the orig.nofix surfaces from the recon-all.log file
+    # see Rosen et al. (2017), https://www.biorxiv.org/content/early/2017/10/01/125161
+    def extract_euler(logfile):
+        with open(logfile) as fi:
+            logtext = fi.read()
+        p = re.compile(r"orig.nofix lheno =\s+(-?\d+), rheno =\s+(-?\d+)")
+        results = p.findall(logtext)
+        if len(results) != 1:
+            raise Exception("Euler number could not be extracted from {}".format(logfile))
+        lh_euler, rh_euler = results[0]
+        return int(lh_euler), int(rh_euler)
+
+
+
+    euler_out_file = os.path.join(table_dir, "euler.tsv")
+    print("Writing euler tables to %s." % euler_out_file)
+
+    # get freesurfer subjects
+    os.chdir(output_dir)
+    subjects = []
+    for s in subjects_to_analyze:
+        subjects += glob("sub-{}*".format(s))
+    # remove long subjects as they don't have orig.nofix surfaces,
+    # therefore no euler numbers
+    subjects = list(filter(lambda s: ".long.sub-" not in s, subjects))
+    if len(subjects) > 0:
+        df = pd.DataFrame([], columns=["subject", "lh_euler", "rh_euler"])
+        for subject in subjects:
+            logfile = os.path.join(output_dir, subject, "scripts/recon-all.log")
+            lh_euler, rh_euler = extract_euler(logfile)
+            df_subject = pd.DataFrame({"subject": [subject],
+                                        "lh_euler": [lh_euler],
+                                        "rh_euler": [rh_euler]},
+                                        columns=["subject", "lh_euler", "rh_euler"])
+            df = df.append(df_subject)
+        df["mean_euler_bh"] = df[["lh_euler", "rh_euler"]].mean(1)
+        df.sort_values("subject", inplace=True)
+        df.to_csv(euler_out_file, sep="\t", index=False)
+    else:
+        print("\nNo subjects included in the analysis. Skipping group2 level euler number table.")