pyBabyAFQ #524

bloomdt-uw · 2020-10-08T06:31:14Z

Based on Mareike's recommendation, pausing work on dHCP dataset and generally abandoning approach to use FA. Instead will focus on replicating MATLAB AFQ results with Kalanit longitudinal data.

Next Steps:

Kalanit longitudinal data

Processed using MATLAB BabyAFQ

Goal is to reproduce results in:

Gold standard
MATLAB AFQ
use existing tractography data and feed into pyAFQ segmentation
- convert dataset to BIDS compliant
- may need to use PyAFQ scripts to convert MATLAB AFQ files
- verify can convert mif files to nii
- possible to use two step process and segment using CLI?
if successful, then consider tractography using ACT within pyAFQ

dHCP data

Single Subject data

if Kalanit pipeline sucessful, continue with dHCP

Revisit BIDS conversion

using default dHCP dataset_description.json generates BIDS validation error.

dHCP uses BIDS 1.0.2 (dHCP Data Release Notes)
- converter from BIDS 1.0.2 to BIDS 1.4.0
- add BIDSVersion to original file
- document and continue to use workaround: replace "PipelineDescription.Name" with `"PipelineDescription": { "Name"
Error processing sub-CC#####XX##_ses-#####_desc-preproc_dwi.json file when present. Currently just excluded file
- investigate error
when making BIDS compliant pyAFQ required data in AFQ_data derivatives folder
- add to documentation
pyAFQ does not read nii files must be nii.gz, see Make s3fs_nifti_read compatible with non-gzipped files #151, therefore needed to gzip nii files.
Convert script to Jupyter notebook with CloudKnot (see [ENH] Add cloudknot example #533)
if pyAFQ tractography doesn't work, consider two-step pipeline: use qsiprep/mrtrix for tractography then pyAFQ segmentation
add segmentation end point filtering
dependency: custom ROI monkey patch can be removed once Make bundle input more flexible #519 merged

debugging FA tractography

However if decide to revisit dHCP FA:

Parameters explored thus far (all without endpoint filtering):

default
b0 1000
b0 1000 with FA seed threshold 0.1
b0 1000 with ROIMask
b0 1000 with BrainMask and ROIMask
b0 1000 with BrainMask, ROIMask, and FA stop threshold 0.1

In all cases some bundles were entirely missing and tractography was highly irregular.

arokem · 2021-03-10T19:25:10Z

examples/plot_baby_afq.py

+    "c405e0dbd9a4091c77b3d1ad200229b4", "ec0aeccc6661d2ee5ed79259383cdcee",
+    "2802cd227b550f6e85df0fec1d515c29", "385addb999dc6d76957d2a35c4ee74bb",
+    "b79f01829bd95682faaf545c72b1d52c", "b79f01829bd95682faaf545c72b1d52c",
+    "e49ba370edca96734d9376f551d413db", "f59e9e69e06325198f70047cd63c3bdc",
+    "ae3bd2931f95adae0280a8f75cd3ca9b", "c409a0036b8c2dd4d03d11fbc6bfbdcd",
+    "c2597a474ea5ec9e3126c35fd238f6b2", "67af59c934147c9f9ff6e0b76c4cc6eb",
+    "72d0bbc0b6162e9291fdc450c103a1f0", "51f5041be63ad0ac10d1ac09e3bf1a8e",
+    "6200f5cdc1402dce46dedd505468b147", "83cb5bf6b9b1eda63c8643871e84a6d4",
+    "2a5d8d309b1256d6e48958e112201c2c", "ba24d0915fdff215a403480d0a7745c9",
+    "1001e833d1344274d543ffd02a66af80", "03e20c5eebcd4d439c4ffb36d26a10d9",
+    "6200f5cdc1402dce46dedd505468b147", "83cb5bf6b9b1eda63c8643871e84a6d4",
+    "a6ae325cce2dc4bb21b52ee4c6ca844f", "a96a31df8f96ccf1c5f871a18e4a2d72",
+    "65b7378ca85689a5f23b1b84a6ed78f0", "ce0d0ea696ef51c671aa118e11192e2d",
+    "ce4918086ca1e829698576bcf95db62b", "96168d2eff74ec2acf9065f499732526",
+    "6b20ba319d44472ec21d6ece937605bb", "26b1cf6ec8bd365dde42e3efe9beeac2",
+    "0b3ccf06564d973bfcfff9a87e74f8b5", "84f3426033a2b225b0920b2622864375",
+    "5351b3cb7efa9aa8e83e266342809ebe", "4e8a34aaba4e0f22a6149f38620dc39d",
+    "682c08f66e8c2cf9e4e60f5ce308c76c", "9077affd4f3a8a1d6b44678cde4b3bf4",
+    "5adf36f00669cc547d5eb978acf46266", "66a8002688ffdf3943722361da90ec6a",
+    "efb5ae138df92019541861f9aa6a4d57", "757ec61078b2e9f9a073871b3216ff7a",
+    "ff1e238c52a21f8cc5d44ac614d9627f", "cf16dd2767c6ab2d3fceb2890f6c3e41",
+    "6016621e244b60b9c69fd44b055e4a03", "fd495a2c94b6b13bfb4cd63e293d3fc0",
+    "bf81a23d80f55e5f1eb0c16717193105",
+    "6f8bf8f70216788d14d9a49a3c664b16",
+    "19df0297d6a2ac21da5e432645d63174",
+]
+
+pediatric_remote_fnames = [
+    "24880625", "24880628", "24880631", "24880634", "24880637", "24880640",
+    "24880643", "24880646", "24880649", "24880652", "24880655", "24880661",
+    "24880664", "24880667", "24880670", "24880673",
+    "24880676", "24880679",
+    "24880685", "24880688",
+    "24880691", "24880694", "24880697", "24880700",
+    "24880703", "24880706", "24880712", "24880715", "24880718", "24880721",
+    "24880724", "24880727", "24880730", "24880733",
+    "24880736", "24880748",
+    "24880739", "24880742",
+    "24880754", "24880757", "24880760", "24880763",
+    "24880769", "24880772", "24880775", "24880778",
+    "24880781", "24880787", "24880790", "24880793", "24880796", "24880802",
+    "24880805", "24880808",
+    "24880616",
+    "24880613",
+    "24986396"
+]
+
+fetch_pediatric_templates = _make_fetcher("fetch_pediatric_templates",
+                                          op.join(afq_home,
+                                                  'pediatric_templates'),
+                                          baseurl, pediatric_remote_fnames,
+                                          pediatric_fnames,
+                                          md5_list=pediatric_md5_hashes,
+                                          doc="Download pediatric templates")
+
+
+# TODO coding: this duplicates read_xxx_template pattern
+# TODO coding: unit testing
+def read_pediatric_templates(resample_to=False):
+    """Load pediatric AFQ templates from file
+
+    Returns
+    -------
+    dict with: keys: names of template ROIs and values: nibabel Nifti1Image
+    objects from each of the ROI nifti files.
+    """
+    files, folder = fetch_pediatric_templates()
+
+    print('Loading pediatric templates...')
+    tic = time.perf_counter()
+
+    pediatric_templates = {}
+    for f in files:
+        img = nib.load(op.join(folder, f))
+        if resample_to:
+            if isinstance(resample_to, str):
+                resample_to = nib.load(resample_to)
+            img = nib.Nifti1Image(reg.resample(img.get_fdata(),
+                                               resample_to,
+                                               img.affine,
+                                               resample_to.affine),
+                                  resample_to.affine)
+        pediatric_templates[f.split('.')[0]] = img
+
+    toc = time.perf_counter()
+    print(f'pediatric templates loaded in {toc - tic:0.4f} seconds')
+
+    # For the arcuate (AF/ARC), reuse the SLF ROIs
+    pediatric_templates['ARC_roi1_L'] = pediatric_templates['SLF_roi1_L']
+    pediatric_templates['ARC_roi1_R'] = pediatric_templates['SLF_roi1_R']
+    pediatric_templates['ARC_roi2_L'] = pediatric_templates['SLFt_roi2_L']
+    pediatric_templates['ARC_roi2_R'] = pediatric_templates['SLFt_roi2_R']
+    pediatric_templates['ARC_roi3_L'] = pediatric_templates['SLFt_roi3_L']
+    pediatric_templates['ARC_roi3_R'] = pediatric_templates['SLFt_roi3_R']
+
+    # For the middle longitudinal fasciculus (MdLF) reuse ILF ROI
+    pediatric_templates['MdLF_roi2_L'] = pediatric_templates['ILF_roi1_L']
+    pediatric_templates['MdLF_roi2_R'] = pediatric_templates['ILF_roi1_R']
+
+    return pediatric_templates


All of this stuff should move into the data module.

arokem · 2021-03-10T19:25:32Z

examples/plot_baby_afq.py

+                bids_path=dhcp_home,
+                dmriprep='dHCP neonatal dMRI pipeline',
+                # bids_filters={"suffix": "dwi"}, # default
+                # custom_tractography_bids_filters=None, # default


We'd like to use a BIDS filter to find already-tracked data here.

Both for the pilot data (from KGS lab) and eventually also for dHCP data (where tracking will be done with qsiprep/mrtrix).

arokem · 2021-03-10T19:27:46Z

examples/plot_baby_afq.py

+def make_pediatric_bundle_dict(bundle_names=pediatric_bundle_names,
+                               seg_algo="afq",
+                               resample_to=False):
+    # pediatric probability maps
+    prob_map_order = ["ATR_L", "ATR_R", "CST_L", "CST_R", "CGC_L", "CGC_R",
+                      "HCC_L", "HCC_R", "FP", "FA", "IFO_L", "IFO_R", "ILF_L",
+                      "ILF_R", "SLF_L", "SLF_R", "UNC_L", "UNC_R",
+                      "ARC_L", "ARC_R", "MdLF_L", "MdLF_R"]
+
+    prob_maps = pediatric_templates['UNCNeo_JHU_tracts_prob-for-babyAFQ']
+    prob_map_data = prob_maps.get_fdata()
+
+    # pediatric bundle dict
+    pediatric_bundles = {}
+
+    # each bundles gets a digit identifier (to be stored in the tractogram)
+    uid = 1
+
+    for name in pediatric_bundle_names:
+        # ROIs that cross the mid-line
+        if name in ["FA", "FP"]:
+            pediatric_bundles[name] = {
+                'ROIs': [pediatric_templates[name + "_L"],
+                         pediatric_templates[name + "_R"],
+                         pediatric_templates["mid-saggital"]],
+                'rules': [True, True, True],
+                'cross_midline': True,
+                'prob_map': prob_map_data[...,
+                                          prob_map_order.index(name)],
+                'uid': uid}
+            uid += 1
+        # SLF is a special case, because it has an exclusion ROI:
+        elif name == "SLF":
+            for hemi in ['_R', '_L']:
+                pediatric_bundles[name + hemi] = {
+                    'ROIs': [pediatric_templates[name + '_roi1' + hemi],
+                             pediatric_templates[name + '_roi2' + hemi],
+                             pediatric_templates["SLFt_roi2" + hemi]],
+                    'rules': [True, True, False],
+                    'cross_midline': False,
+                    'prob_map': prob_map_data[...,
+                                              prob_map_order.index(name + hemi)],
+                    'uid': uid}
+                uid += 1
+        # Third ROI for curvy tracts
+        elif name in ["ARC", "ATR", "CGC", "IFO", "UNC"]:
+            for hemi in ['_R', '_L']:
+                pediatric_bundles[name + hemi] = {
+                    'ROIs': [pediatric_templates[name + '_roi1' + hemi],
+                             pediatric_templates[name + '_roi2' + hemi],
+                             pediatric_templates[name + '_roi3' + hemi]],
+                    'rules': [True, True, True],
+                    'cross_midline': False,
+                    'prob_map': prob_map_data[...,
+                                              prob_map_order.index(name + hemi)],
+                    'uid': uid}
+                uid += 1
+        elif name == "MdLF":
+            for hemi in ['_R', '_L']:
+                pediatric_bundles[name + hemi] = {
+                    'ROIs': [pediatric_templates[name + '_roi1' + hemi],
+                             pediatric_templates[name + '_roi2' + hemi]],
+                    'rules': [True, True],
+                    'cross_midline': False,
+                    # reuse probability map from ILF
+                    'prob_map': prob_map_data[...,
+                                              prob_map_order.index("ILF" + hemi)],
+                    'uid': uid}
+                uid += 1
+        # Default: two ROIs within hemisphere
+        else:
+            for hemi in ['_R', '_L']:
+                pediatric_bundles[name + hemi] = {
+                    'ROIs': [pediatric_templates[name + '_roi1' + hemi],
+                             pediatric_templates[name + '_roi2' + hemi]],
+                    'rules': [True, True],
+                    'cross_midline': False,
+                    'prob_map': prob_map_data[...,
+                                              prob_map_order.index(name + hemi)],
+                    'uid': uid}
+                uid += 1
+
+    return pediatric_bundles


This probably also moves into data or possibly into baby.data

arokem · 2021-03-10T19:33:47Z

One possible API here is that we create a new module AFQ.baby that includes all the baby-specific stuff.

The question is what do we do about things that are implemented in AFQ.api? Should sub-class AFQ.api.AFQ into AFQ.baby.api.AFQ? We would then need to replace some of the stuff related to registration.

arokem · 2021-03-10T19:39:20Z

examples/plot_baby_afq.py

+                max_bval=1000,  # override
+                reg_template=pediatric_templates['UNCNeo-withCerebellum-for-babyAFQ'],
+                # reg_subject="power_map", # default
+                reg_subject="b0",  # override


Both of these (reg_template and reg_subject) should be skull-stripped using a brain-mask (for pilot data, we have that, for dHCP data it also comes with the anatomical pipeline derivatives). We should over-ride the pyAFQ default brain-mask (which is generated using the median_otsu algorithm) with the brain-mask that is provided with the respective datasets.

arokem · 2021-03-10T19:41:49Z

examples/plot_baby_afq.py

+                reg_subject="b0",  # override
+                # brain_mask=B0Mask(), # default
+                brain_mask=MaskFile("brainmask",
+                                    {"scope": "dHCP neonatal dMRI pipeline"}),


Yep - this is it!

Check if this is a binary mask or a continuous probability mask, and then threshold if needed (there is a class for that!).

Just so I understand what is going on: Which registration algorithm are we using? Are we applying a brain mask in the step above or are we just checking that it is not a probability mask?

reg_algo=None by default, which should use Symmetric Diffeomorphic Registration (SyN) Algorithm for registration.

In this step we are configuring pyAFQ and nothing is actually run until the myafq.export_all() command. According to the docs the brain_mask is applied before registration to a template.

arokem · 2021-03-10T19:49:09Z

examples/plot_baby_afq.py

+                # min_bval=None, # default
+                min_bval=1000,  # override
+                # max_bval=None, # default
+                max_bval=1000,  # override


We might consider adding all the shells here. Proceed with caution.

Eventually, we might have some more denoising coming (through #625), so we could do DKI with multi-shell, and add some more scalars to the profiles (AWF, MK).

arokem · 2021-03-10T19:53:16Z

if successful, then consider tractography using ACT within pyAFQ

Let's punt on this one for now. We'll use the tractography that is already there for the pilot data and do qsiprep/mrtrix with the dHCP data.

arokem · 2021-03-10T19:57:43Z

add segmentation end point filtering

Might not need this for this application. The MATLAB babyAFQ does not use endpoints for filtering.

arokem · 2021-03-10T20:00:56Z

add segmentation end point filtering

Furthermore, the template does have a version of the AAL atlas, so if need be, we can come back to that (with the caveat that we'd need to do some massaging of the template to add some volumes that we have in the template we use with grownups).

We'll also need the ALS atlas for the definition of pAF (which, in this PR, is currently being skipped). In this case, a particular part of the AAL atlas is extracted to serve as a waypoint ROI - this is already part of the templates.

pep8speaks · 2021-03-11T00:38:39Z

Hello @bloomdt-uw! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

In the file AFQ/api.py:

Line 325:80: E501 line too long (89 > 79 characters)
Line 338:1: W293 blank line contains whitespace
Line 344:80: E501 line too long (87 > 79 characters)
Line 346:80: E501 line too long (84 > 79 characters)
Line 351:14: E114 indentation is not a multiple of four (comment)
Line 351:14: E116 unexpected indentation (comment)
Line 352:80: E501 line too long (83 > 79 characters)
Line 353:80: E501 line too long (86 > 79 characters)
Line 356:1: W293 blank line contains whitespace
Line 357:80: E501 line too long (81 > 79 characters)
Line 363:80: E501 line too long (83 > 79 characters)
Line 371:33: E128 continuation line under-indented for visual indent
Line 372:33: E128 continuation line under-indented for visual indent
Line 376:49: E128 continuation line under-indented for visual indent
Line 383:80: E501 line too long (80 > 79 characters)
Line 384:37: E128 continuation line under-indented for visual indent
Line 385:37: E128 continuation line under-indented for visual indent
Line 389:53: E128 continuation line under-indented for visual indent
Line 389:80: E501 line too long (87 > 79 characters)
Line 396:80: E501 line too long (80 > 79 characters)
Line 397:37: E128 continuation line under-indented for visual indent
Line 398:37: E128 continuation line under-indented for visual indent
Line 398:80: E501 line too long (80 > 79 characters)
Line 402:53: E128 continuation line under-indented for visual indent
Line 402:80: E501 line too long (87 > 79 characters)
Line 408:80: E501 line too long (80 > 79 characters)
Line 409:37: E128 continuation line under-indented for visual indent
Line 409:80: E501 line too long (80 > 79 characters)
Line 414:53: E128 continuation line under-indented for visual indent
Line 414:80: E501 line too long (88 > 79 characters)
Line 421:80: E501 line too long (80 > 79 characters)
Line 422:37: E128 continuation line under-indented for visual indent
Line 422:80: E501 line too long (80 > 79 characters)
Line 426:53: E128 continuation line under-indented for visual indent
Line 426:80: E501 line too long (87 > 79 characters)
Line 1075:1: W293 blank line contains whitespace

In the file AFQ/data.py:

Line 205:80: E501 line too long (138 > 79 characters)
Line 300:36: W291 trailing whitespace
Line 301:1: W293 blank line contains whitespace
Line 306:11: W291 trailing whitespace
Line 307:43: W291 trailing whitespace

Comment last updated at 2021-08-03 22:23:15 UTC

arokem

I think that this is really close to being ready to merge! Just a few small comments here.

arokem · 2021-08-03T21:14:38Z

AFQ/api.py

@@ -945,6 +1070,9 @@ def _get_best_scalar(self):
        return self.scalars[0]

    def get_reg_template(self):
+        if isinstance(self.reg_template, nib.Nifti1Image):


@36000 : what do you think? Does this change make sense? Or should this go somewhere else?

@36000 : do you have any thoughts here? I would like to merge this PR, but need some input from you for this particular line.

arokem · 2021-08-03T21:15:50Z

examples/plot_baby_afq.py

+
+The following is an example of tractometry for pediatric bundles.
+
+TODO cite @mareike


https://www.biorxiv.org/content/10.1101/2021.03.29.437583v1

arokem · 2021-08-03T21:43:52Z

docs/source/conf.py

@@ -209,6 +209,7 @@
     'examples_dirs': '../../examples',
     # path where to save gallery generated examples
     'gallery_dirs': 'auto_examples',
+     'ignore_pattern': 'plot_baby_afq.py',


I wonder why this is not in the exclude_patterns on line 82 of this file (and what the difference is between these), but maybe that's not a specific issue with this PR

I don't recall exactly why I made this change, it was part of the original commit, so it could be OBE

arokem · 2021-08-03T21:45:54Z

examples/plot_baby_afq.py

+# -----------------
+# **Diffusion dataset**
+#
+# .. note::


We don't actually use any dataset here, so I would say that this is all somewhat theoretical, and link to the repo that contains the code for the dHCP data analysis, to say that's what we actually do.

arokem · 2021-08-04T22:37:22Z

I think this is good to go! @36000 : would you mind taking a quick look to confirm?

arokem · 2021-09-26T20:03:06Z

OK - I believe that the change alluded to here is not harmful, so I am going to go ahead and merge this, so that we can avoid conflicts down the line.

bloomdt-uw force-pushed the enh-487-BabyAFQ branch 3 times, most recently from 07e4908 to 27eea49 Compare October 10, 2020 00:18

arokem reviewed Mar 10, 2021

View reviewed changes

bloomdt-uw force-pushed the enh-487-BabyAFQ branch from 341e259 to 69cf219 Compare March 11, 2021 02:08

bloomdt-uw force-pushed the enh-487-BabyAFQ branch from 69cf219 to 18ac9dd Compare March 21, 2021 18:52

bloomdt-uw force-pushed the enh-487-BabyAFQ branch from 18ac9dd to 638de57 Compare June 3, 2021 16:37

bloomdt-uw added 3 commits August 2, 2021 11:14

Initial pyBabyAFQ

e65e9bf

disable sphinx-build for babyafq and adjusting afq pipeline parameters

6e843c1

reuse ILF probability map for MdLF

83a3173

bloomdt-uw force-pushed the enh-487-BabyAFQ branch from 638de57 to 83a3173 Compare August 2, 2021 18:15

bloomdt-uw added 2 commits August 3, 2021 13:45

adding pedatric templates and bundle specifications to pyAFQ

7786030

removing unnecessary condition

4ca9367

arokem changed the title ~~WIP pyBabyAFQ~~ pyBabyAFQ Aug 3, 2021

arokem marked this pull request as ready for review August 3, 2021 21:40

arokem reviewed Aug 3, 2021

View reviewed changes

updating example documentation

0cf2310

arokem merged commit 7961e2b into yeatmanlab:master Sep 26, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

pyBabyAFQ #524

pyBabyAFQ #524

bloomdt-uw commented Oct 8, 2020 •

edited

arokem Mar 10, 2021

arokem Mar 10, 2021

arokem Mar 10, 2021

arokem commented Mar 10, 2021

arokem Mar 10, 2021

arokem Mar 10, 2021

grotheer Mar 15, 2021

bloomdt-uw Mar 16, 2021

arokem Mar 10, 2021

arokem commented Mar 10, 2021

arokem commented Mar 10, 2021

arokem commented Mar 10, 2021

pep8speaks commented Mar 11, 2021 •

edited

arokem left a comment

arokem Aug 3, 2021

arokem Aug 12, 2021

arokem Aug 3, 2021

arokem Aug 3, 2021

bloomdt-uw Aug 3, 2021

arokem Aug 3, 2021

arokem commented Aug 4, 2021

arokem commented Sep 26, 2021


		The following is an example of tractometry for pediatric bundles.

		TODO cite @mareike

pyBabyAFQ #524

pyBabyAFQ #524

Conversation

bloomdt-uw commented Oct 8, 2020 • edited

Next Steps:

debugging FA tractography

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Mar 10, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Mar 10, 2021

arokem commented Mar 10, 2021

arokem commented Mar 10, 2021

pep8speaks commented Mar 11, 2021 • edited

Comment last updated at 2021-08-03 22:23:15 UTC

arokem left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Aug 4, 2021

arokem commented Sep 26, 2021

bloomdt-uw commented Oct 8, 2020 •

edited

pep8speaks commented Mar 11, 2021 •

edited