Getting started with pyAFQ - GroupAFQ#

There are two ways to :doc:use pyAFQ </tutorials/index>: through the command line interface, and by writing Python code. This tutorial will walk you through the basics of the latter, using pyAFQ’s Python Application Programming Interface (API).

import os.path as op

import matplotlib.pyplot as plt
import nibabel as nib
import plotly
import pandas as pd

from AFQ.api.group import GroupAFQ
import AFQ.data.fetch as afd
import AFQ.viz.altair as ava
import AFQ.definitions.image as afm
2026-05-19 01:03:38,878	INFO util.py:154 -- Missing packages: ['ipywidgets']. Run `pip install -U ipywidgets`, then restart the notebook server for rich notebook output.

Example data#

pyAFQ can be called using GroupAFQ to handle data organized in a BIDS compliant directory. If this is not the case, refer to the Participant AFQ example. To get users started with this tutorial, we will download some example data and organize it in a BIDS compliant way (for more details on how BIDS is used in pyAFQ, refer to :doc:plot_006_bids_layout).

The following call downloads a a single subject’s data from the Healthy Brain Network Processed Open Diffusion Derivatives dataset (HBN-POD2) [1], [2] and organizes it in BIDS in the user’s home directory under::

~/AFQ_data/HBN/

The data is also placed in a derivatives directory, signifying that it has already undergone the required preprocessing necessary for pyAFQ to run.

The clear_previous_afq is used to remove any previous runs of the afq object stored in the ~/AFQ_data/HBN/ BIDS directory. Set it to None if you want to use the results of previous runs.

bids_path = afd.fetch_hbn_preproc(
    ["NDARAA948VFH"],
    clear_previous_afq="all")[1]

Set tractography parameters (optional)#

We make create a tracking_params variable, which we will pass to the GroupAFQ object which specifies that we want 100,000 seeds randomly distributed in the white matter. We only do this to make this example faster and consume less space; normally, we use more seeds.

tracking_params = dict(n_seeds=1e5,
                       random_seeds=True,
                       rng_seed=2025,
                       trx=True)

Define PVE images (optional)#

To improve segmentation and tractography results, we can provide partial volume estimate (PVE) images for the cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM). Here, we define these images using the AFQ.definitions.image.PVEImages class, which takes as input three AFQ.definitions.image.ImageFile objects, one for each tissue type. One can also provide a single PVE image with all three tissue types using the AFQ.definitions.image.PVEImage class. Finally, by default, if no PVE images are provided, pyAFQ will use SynthSeg2 to compute these images.

pve = afm.PVEImages(
    afm.ImageFile(
        suffix="probseg", filters={"scope": "qsiprep", "label": "CSF"}),
    afm.ImageFile(
        suffix="probseg", filters={"scope": "qsiprep", "label": "GM"}),
    afm.ImageFile(
        suffix="probseg", filters={"scope": "qsiprep", "label": "WM"}))

Brain Mask Definition (optional)#

By default, pyAFQ will compute a brain mask from the T1. However, this requires onnxruntime to be installed. If you do not have onnxruntime installed, or if you want to use a different brain mask, you can specify it here.

brain_mask_definition = afm.ImageFile(
    suffix="mask", filters={"desc": "brain", "scope": "qsiprep"})

Initialize a GroupAFQ object:#

Creates a GroupAFQ object, that encapsulates tractometry. This object can be used to manage the entire :doc:/explanations/index, including:

  • Tractography

  • Registration

  • Segmentation

  • Cleaning

  • Profiling

  • Visualization

This will also create an output folder for the corresponding AFQ derivatives in the AFQ data directory: AFQ_data/HBN/derivatives/afq/

To initialize this object we will pass in the path location to our BIDS compliant data, the name of the preprocessing pipeline we want to use, the name of the t1 preprocessing pipeline we want to use (in this case, its the same, qsiprep [3]), the participant labels we want to process (in this case, just a single subject), the PVE images we defined above, and the tracking parameters we defined above.

myafq = GroupAFQ(
    bids_path=op.join(afd.afq_home, 'HBN'),
    dwi_preproc_pipeline='qsiprep',
    t1_preproc_pipeline='qsiprep',
    participant_labels=['NDARAA948VFH'],
    pve=pve,
    brain_mask_definition=brain_mask_definition,
    tracking_params=tracking_params)
INFO:bidsschematools:No schema path specified, defaulting to the bundled schema, `/opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/bidsschematools/data/schema.json`.
INFO:AFQ:Using the following files for subject NDARAA948VFH and session HBNsiteRU:
INFO:AFQ:  DWI: /home/runner/AFQ_data/HBN/derivatives/qsiprep/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.nii.gz
INFO:AFQ:  BVAL: /home/runner/AFQ_data/HBN/derivatives/qsiprep/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.bval
INFO:AFQ:  BVEC: /home/runner/AFQ_data/HBN/derivatives/qsiprep/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.bvec
INFO:AFQ:  T1: /home/runner/AFQ_data/HBN/derivatives/qsiprep/sub-NDARAA948VFH/anat/sub-NDARAA948VFH_desc-preproc_T1w.nii.gz

Calculating DKI FA (Diffusion Kurtosis Imaging Fractional Anisotropy)#

The GroupAFQ object has a method called export, which allows the user to calculate various derived quantities from the data.

For example, FA can be computed using the DKI model, by explicitly calling myafq.export("dki_fa"). This triggers the computation of DKI parameters for all subjects in the dataset, and stores the results in the AFQ derivatives directory. In addition, it calculates the FA from these parameters and stores it in a different file in the same directory.

Note

The AFQ API computes quantities lazily. This means that DKI parameters are not computed until they are required. This means that the first line below is the one that requires time.

The result of the call to export is a dictionary, with the subject IDs as keys, and the filenames of the corresponding files as values. This means that to extract the filename corresponding to the FA of the first subject, we can do:

FA_fname = myafq.export("dki_fa", collapse=False)["NDARAA948VFH"]["HBNsiteRU"]

# We will then use `nibabel` to load the deriviative file and retrieve the
# data array.

FA_img = nib.load(FA_fname)
FA = FA_img.get_fdata()
INFO:AFQ:Calculating _b0ref.nii.gz...
INFO:AFQ:_b0ref.nii.gz completed. Saving to /home/runner/AFQ_data/HBN/derivatives/afq/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_b0ref.nii.gz
INFO:AFQ:Calculating _desc-brain_mask.nii.gz...
INFO:AFQ:Calculating _desc-T1w_mask.nii.gz...
INFO:AFQ:_desc-T1w_mask.nii.gz completed. Saving to /home/runner/AFQ_data/HBN/derivatives/afq/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_desc-T1w_mask.nii.gz
INFO:AFQ:_desc-brain_mask.nii.gz completed. Saving to /home/runner/AFQ_data/HBN/derivatives/afq/sub-NDARAA948VFH/ses-HBNsiteRU/dwi/sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_desc-brain_mask.nii.gz
INFO:AFQ:Calculating _model-kurtosis_param-diffusivity_dwimap.nii.gz, _model-kurtosis_param-s0_dwimap.nii.gz...
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[7], line 1
----> 1 FA_fname = myafq.export("dki_fa", collapse=False)["NDARAA948VFH"]["HBNsiteRU"]
      2 
      3 # We will then use `nibabel` to load the deriviative file and retrieve the
      4 # data array.

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/api/group.py:647, in GroupAFQ.export(self, attr_name, collapse)
    645         to_calc_list.append((subject, session))
    646     else:
--> 647         results[subject][session] = val_from_plan(plans_dict, attr_name)
    649 # if some need to be calculated, do those in parallel
    650 if to_calc_list:

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/api/utils.py:146, in val_from_plan(plan, attr_name)
    144 def val_from_plan(plan, attr_name):
    145     try:
--> 146         return plan[attr_name]
    147     except Exception as err:
    148         current_err = err

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/pcollections/_lazy.py:381, in ldict.__getitem__(self, key)
    379 v = pdict.__getitem__(self, key)
    380 if isinstance(v, lazy):
--> 381     v = v()
    382 return v

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/pcollections/_lazy.py:143, in lazy.__call__(self)
    141     val = self.value
    142 else:
--> 143     val = part()
    144     # We've successfully calculated the value; set the members
    145     # appropriately.
    146     object.__setattr__(self, 'value', val)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:855, in plan._source_lookup(inputtup, calctup, src)
    853     (cidx, oidx) = src
    854     lazycalc = calctup[cidx]
--> 855     val = lazycalc()[oidx]
    856 else:
    857     val = inputtup[src]()

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/pcollections/_lazy.py:143, in lazy.__call__(self)
    141     val = self.value
    142 else:
--> 143     val = part()
    144     # We've successfully calculated the value; set the members
    145     # appropriately.
    146     object.__setattr__(self, 'value', val)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:884, in plan._make_calctup.<locals>.<lambda>(c, args)
    881 # We take advantage of Python's weak closures here:
    882 calctup = ()
    883 calctup = tuple(
--> 884     lazy(lambda c,args: f(inputtup, calctup, c, args), c, args)
    885     for (c,args) in zip(calcdata.calcs, calcdata.args))
    886 return calctup

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:868, in plan._call_calc(inputtup, calctup, c, args)
    866 kwargs = {}
    867 c = to_calc(c)
--> 868 for (p,arg) in zip(c.signature.parameters.values(), argvals):
    869     if p.kind == p.POSITIONAL_ONLY:
    870         args.append[arg]

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:855, in plan._source_lookup(inputtup, calctup, src)
    853     (cidx, oidx) = src
    854     lazycalc = calctup[cidx]
--> 855     val = lazycalc()[oidx]
    856 else:
    857     val = inputtup[src]()

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/pcollections/_lazy.py:143, in lazy.__call__(self)
    141     val = self.value
    142 else:
--> 143     val = part()
    144     # We've successfully calculated the value; set the members
    145     # appropriately.
    146     object.__setattr__(self, 'value', val)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:884, in plan._make_calctup.<locals>.<lambda>(c, args)
    881 # We take advantage of Python's weak closures here:
    882 calctup = ()
    883 calctup = tuple(
--> 884     lazy(lambda c,args: f(inputtup, calctup, c, args), c, args)
    885     for (c,args) in zip(calcdata.calcs, calcdata.args))
    886 return calctup

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:873, in plan._call_calc(inputtup, calctup, c, args)
    871     else:
    872         kwargs[p.name] = arg
--> 873 r = c.eager_call(*args, **kwargs)
    874 if is_amap(r):
    875     return tuple(map(r.__getitem__, c.outputs))

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/immlib/workflow/_core.py:448, in calc.eager_call(self, *args, **kwargs)
    432 """Eagerly calls the given calculation using the arguments.
    433 
    434 ``c.eager_call(...)`` returns the result of calling the calculation
   (...)    444 calc.eager_mapcall, calc.lazy_call, calc.lazy_mapcall
    445 """
    446 # Now we just pass these arguments along (the function itself has been
    447 # given the caching code via decorators already).
--> 448 res = self.function(*args, **kwargs)
    449 # Now interpret the result.
    450 outs = self.outputs

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/tasks/decorators.py:104, in as_file.<locals>._as_file.<locals>.wrapper_as_file(*args, **kwargs)
    101 logger.info(f"Calculating {calculation_name}...")
    103 try:
--> 104     results = func(*args, **kwargs)
    106     if len(output_specs) == 1:
    107         results = [results]

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/tasks/decorators.py:228, in as_img.<locals>.wrapper_as_img(*args, **kwargs)
    225     dwi_affine = data_imap["dwi_affine"]
    227 start_time = time()
--> 228 results = func(*args, **kwargs)
    229 elapsed = time() - start_time
    231 is_single_output = isinstance(results[0], np.ndarray)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/tasks/data.py:318, in dki_params(brain_mask, gtab, data, citations)
    310     raise ValueError(
    311         (
    312             "The DKI model requires at least 2 non-zero b-values, "
   (...)    315         )
    316     )
    317 mask = nib.load(brain_mask).get_fdata()
--> 318 dkf = dki_fit_model(gtab, data, mask=mask, return_S0_hat=True)
    319 meta = dict(
    320     Description=(
    321         "Diffusion Coefficient, encoded as a kurtosis tensor representation"
   (...)    330     ),
    331 )
    333 meta_s0 = dict(
    334     Description="Estimated signal intensity with no diffusion weighting, ie. S0",
    335     Model=dict(
   (...)    341     ),
    342 )

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/models/dki.py:16, in _fit(gtab, data, mask, return_S0_hat)
     14 def _fit(gtab, data, mask=None, return_S0_hat=False):
     15     dkimodel = dki.DiffusionKurtosisModel(gtab, return_S0_hat=return_S0_hat)
---> 16     return dkimodel.fit(
     17         data,
     18         mask=mask,
     19     )

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:201, in warning_for_keywords.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    194 # Check if the current version is within the warning range
    195 if (
    196     version.parse(from_version)
    197     <= version.parse(current_version)
    198     <= version.parse(until_version)
    199 ):
    200     # Convert positional to keyword arguments and issue a warning
--> 201     return convert_positional_to_keyword(func, args, kwargs)
    203 # If the version is greater than the until_version,
    204 # pass the arguments as they are
    205 elif version.parse(current_version) > version.parse(until_version):

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:192, in warning_for_keywords.<locals>.decorator.<locals>.wrapper.<locals>.convert_positional_to_keyword(func, args, kwargs)
    182         warnings.warn(
    183             f"Pass {positionally_passed_kwonly_args} as keyword args. "
    184             f"From version {until_version} passing these as positional "
   (...)    187             stacklevel=3,
    188         )
    190     return func(*positional_args, **corrected_kwargs)
--> 192 return func(*args, **kwargs)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/reconst/dki.py:1896, in DiffusionKurtosisModel.fit(self, data, mask)
   1893 data_thres = np.maximum(data, self.min_signal)
   1895 if self.is_multi_method and not self.is_iter_method:
-> 1896     fit_result, extra = self.multi_fit(
   1897         data_thres, mask=mask, weights=self.weights, **self.kwargs
   1898     )
   1899     if extra is not None:
   1900         self.extra = extra

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/reconst/multi_voxel.py:282, in multi_voxel_fit.<locals>.decorator.<locals>.new_fit(self, data, mask, **kwargs)
    279 if weights_is_array:
    280     kwargs["weights"] = weights[ijk]
--> 282 svf = single_voxel_fit(self, data[ijk], **kwargs)
    284 # Not all fit methods return extra, handle this here
    285 if isinstance(svf, tuple):

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:201, in warning_for_keywords.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    194 # Check if the current version is within the warning range
    195 if (
    196     version.parse(from_version)
    197     <= version.parse(current_version)
    198     <= version.parse(until_version)
    199 ):
    200     # Convert positional to keyword arguments and issue a warning
--> 201     return convert_positional_to_keyword(func, args, kwargs)
    203 # If the version is greater than the until_version,
    204 # pass the arguments as they are
    205 elif version.parse(current_version) > version.parse(until_version):

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:192, in warning_for_keywords.<locals>.decorator.<locals>.wrapper.<locals>.convert_positional_to_keyword(func, args, kwargs)
    182         warnings.warn(
    183             f"Pass {positionally_passed_kwonly_args} as keyword args. "
    184             f"From version {until_version} passing these as positional "
   (...)    187             stacklevel=3,
    188         )
    190     return func(*positional_args, **corrected_kwargs)
--> 192 return func(*args, **kwargs)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/reconst/dki.py:1986, in DiffusionKurtosisModel.multi_fit(self, data, mask, **kwargs)
   1976 """Convenience function for fitting multiple voxels."""
   1977 extra_args = (
   1978     {}
   1979     if not self.convexity_constraint
   (...)   1983     }
   1984 )
-> 1986 params, extra = self.fit_method(
   1987     self.design_matrix,
   1988     data,
   1989     self.inverse_design_matrix,
   1990     return_S0_hat=self.return_S0_hat,
   1991     min_diffusivity=self.min_diffusivity,
   1992     **extra_args,
   1993     **kwargs,
   1994 )
   1996 S0_params = None
   1997 if self.return_S0_hat:

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:201, in warning_for_keywords.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    194 # Check if the current version is within the warning range
    195 if (
    196     version.parse(from_version)
    197     <= version.parse(current_version)
    198     <= version.parse(until_version)
    199 ):
    200     # Convert positional to keyword arguments and issue a warning
--> 201     return convert_positional_to_keyword(func, args, kwargs)
    203 # If the version is greater than the until_version,
    204 # pass the arguments as they are
    205 elif version.parse(current_version) > version.parse(until_version):

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:192, in warning_for_keywords.<locals>.decorator.<locals>.wrapper.<locals>.convert_positional_to_keyword(func, args, kwargs)
    182         warnings.warn(
    183             f"Pass {positionally_passed_kwonly_args} as keyword args. "
    184             f"From version {until_version} passing these as positional "
   (...)    187             stacklevel=3,
    188         )
    190     return func(*positional_args, **corrected_kwargs)
--> 192 return func(*args, **kwargs)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/reconst/dki.py:2793, in ls_fit_dki(design_matrix, data, inverse_design_matrix, return_S0_hat, weights, min_diffusivity, return_lower_triangular, return_leverages)
   2790     return result, leverages
   2792 # Write output
-> 2793 dki_params = params_to_dki_params(result, min_diffusivity=min_diffusivity)
   2795 if return_S0_hat:
   2796     return (dki_params[..., 0:-1], dki_params[..., -1]), leverages

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:201, in warning_for_keywords.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    194 # Check if the current version is within the warning range
    195 if (
    196     version.parse(from_version)
    197     <= version.parse(current_version)
    198     <= version.parse(until_version)
    199 ):
    200     # Convert positional to keyword arguments and issue a warning
--> 201     return convert_positional_to_keyword(func, args, kwargs)
    203 # If the version is greater than the until_version,
    204 # pass the arguments as they are
    205 elif version.parse(current_version) > version.parse(until_version):

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:192, in warning_for_keywords.<locals>.decorator.<locals>.wrapper.<locals>.convert_positional_to_keyword(func, args, kwargs)
    182         warnings.warn(
    183             f"Pass {positionally_passed_kwonly_args} as keyword args. "
    184             f"From version {until_version} passing these as positional "
   (...)    187             stacklevel=3,
    188         )
    190     return func(*positional_args, **corrected_kwargs)
--> 192 return func(*args, **kwargs)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/reconst/dki.py:2677, in params_to_dki_params(result, min_diffusivity)
   2675 # Extracting the diffusion tensor parameters from solution
   2676 DT_elements = result[:6]
-> 2677 evals, evecs = decompose_tensor(
   2678     from_lower_triangular(DT_elements), min_diffusivity=min_diffusivity
   2679 )
   2681 # Extracting kurtosis tensor parameters from solution
   2682 MD_square = evals.mean(0) ** 2

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:196, in warning_for_keywords.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    192     return func(*args, **kwargs)
    194 # Check if the current version is within the warning range
    195 if (
--> 196     version.parse(from_version)
    197     <= version.parse(current_version)
    198     <= version.parse(until_version)
    199 ):
    200     # Convert positional to keyword arguments and issue a warning
    201     return convert_positional_to_keyword(func, args, kwargs)
    203 # If the version is greater than the until_version,
    204 # pass the arguments as they are

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/packaging/version.py:122, in parse(version)
    111 def parse(version: str) -> Version:
    112     """Parse the given version string.
    113 
    114     This is identical to the :class:`Version` constructor.
   (...)    120     :raises InvalidVersion: When the version string is not a valid version.
    121     """
--> 122     return Version(version)

File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/packaging/version.py:414, in Version.__init__(self, version)
    404 def __init__(self, version: str) -> None:
    405     """Initialize a Version object.
    406 
    407     :param version:
   (...)    412         exception will be raised.
    413     """
--> 414     if _SIMPLE_VERSION_INDICATORS.issuperset(version):
    415         try:
    416             self._release = tuple(map(int, version.split(".")))

KeyboardInterrupt: 

Visualize the result with Matplotlib#

At this point FA is an array, and we can use standard Python tools to visualize it or perform additional computations with it.

In this case we are going to take an axial slice halfway through the FA data array and plot using a sequential color map.

Note

The data array is structured as a xyz coordinate system.

fig, ax = plt.subplots(1)
ax.matshow(FA[:, :, FA.shape[-1] // 2], cmap='viridis')
ax.axis("off")

Recognizing the bundles and calculating tract profiles:#

Typically, users of pyAFQ are interested in calculating not only an overall map of the FA, but also the major white matter pathways (or bundles) and tract profiles of tissue properties along their length. To trigger the pyAFQ pipeline that calculates the profiles, users can call the export('profiles') method:

Note

Running the code below triggers the full pipeline of operations leading to the computation of the tract profiles. Therefore, it takes a little while to run (about 40 minutes, typically).

myafq.export('profiles')

Visualizing the bundles and calculating tract profiles:#

The pyAFQ API provides several ways to visualize bundles and profiles.

First, we will run a function that exports an html file that contains an interactive visualization of the bundles that are segmented.

Note

By default we resample a 100 points within a bundle, however to reduce processing time we will only resample 50 points.

Once it is done running, it should pop a browser window open and let you interact with the bundles.

Note

You can hide or show a bundle by clicking the legend, or select a single bundle by double clicking the legend. The interactive visualization will also all you to pan, zoom, and rotate.

bundle_html = myafq.export("all_bundles_figure", collapse=False)
plotly.io.show(bundle_html["NDARAA948VFH"]["HBNsiteRU"][0])

We can also visualize the tract profiles in all of the bundles. These plots show both FA (left) and MD (right) laid out anatomically. To make this plot, it is required that you install with pip install pyAFQ[plot] so that you have the necessary dependencies.

fig_files = myafq.export("tract_profile_plots", collapse=False)[
    "NDARAA948VFH"]["HBNsiteRU"]

.. figure:: {{ fig_files[0] }}

We can even use altair to visualize the tract profiles in all of the bundles. We provide a more customizable interface for visualizing the tract profiles using altair. Again, to make this plot, it is required that you install with pip install pyAFQ[plot] so that you have the necessary dependencies.

profiles_df = myafq.combine_profiles()
altair_df = ava.combined_profiles_df_to_altair_df(
    profiles_df,
    tissue_properties=['dki_fa', 'dki_md'])
altair_chart = ava.altair_df_to_chart(altair_df)
altair_chart.display()

Exporting citations#

Finally, we can export the citations for the some of methods used in this analysis. These are not guaranteed to be comprehensive, but they should be a good starting point.

myafq.export("citations")

We can check the number of streamlines per bundle, to make sure every bundle is found with a reasonable amount of streamlines.

bundle_counts = pd.read_csv(
    myafq.export("sl_counts", collapse=False)[
        "NDARAA948VFH"]["HBNsiteRU"], index_col=[0])
for ind in bundle_counts.index:
    if ind == "Total Recognized":
        threshold = 2500
    else:
        threshold = 20
    if bundle_counts["n_streamlines"][ind] < threshold:
        raise ValueError((
            "Small number of streamlines found "
            f"for bundle(s):\n{bundle_counts}"))

Alternative way to initialize GroupAFQ when using QSIPrep data#

As a final note, if you are using QSIPrep preprocessed data, you can also initialize the GroupAFQ object using the from_qsiprep class method. This method will automatically set the appropriate BIDS filters to find the preprocessed DWI data. Additionally, it will find and use the brain masks and PVE images that QSIPrep generates. Outside of BIDS filters, the arguments are the same as those used when initializing the GroupAFQ object directly.

myafq = GroupAFQ.from_qsiprep(
    qsi_dir=op.join(afd.afq_home, 'HBN'),
    participant_labels=['NDARAA948VFH'],
    tracking_params=tracking_params)

References#

.. [1] Alexander LM, Escalera J, Ai L, et al. An open resource for transdiagnostic research in pediatric mental health and learning disorders. Sci Data. 2017;4:170181.

.. [2] Richie-Halford A, Cieslak M, Ai L, et al. An analysis-ready and quality controlled resource for pediatric brain white-matter research. Scientific Data. 2022;9(1):1-27.

.. [3] Cieslak M, Cook PA, He X, et al. QSIPrep: an integrative platform for preprocessing and reconstructing diffusion MRI data. Nat Methods. 2021;18(7):775-778.