Substantially Fewer Steamlines and Tracts using v0.11/v0.12 compared to v0.8.

[jdias001]

The output tractograms I am getting with version 0.12 are very different from those I had in version 0.8. In particular, the number of fibers and tracts identified is substantially less than before and results in very little usable data.

The code I used in v0.8:

from AFQ.data import fetch
from AFQ.api.group import GroupAFQ
import AFQ.definitions.mask as afm
import os
from glob import glob

myafq = api.AFQ(bids_path="C:\\Users\\james\\Documents\\DKI\\pyAFQ\\BIDS\\dmriprep\\",
                scalars=["dti_fa","dti_md","dki_fa","dki_md","dki_awf","dki_mk"],
                tracking_params=dict(
                    odf_model="DTI",
                    ),
                segmentation_params=dict(
                    seg_algo="afq"
                    ),
               )
myafq.export_all()

The output I get is:

However, using code from v0.12:
(Forgive the additional code for looping individual participant data).

import os
from glob import glob
import shutil
import re

from AFQ.api.participant import ParticipantAFQ
from AFQ.definitions.image import ImageFile

BIDSfolder = "C:\\Users\\james\\Documents\\DKI\\pyAFQ\\BIDS\\dmriprep\\"
infolder = "C:\\Users\\james\\Documents\\DKI\\pyAFQ\\BIDS\\derivatives\\pydesigner\\"
outfolder = "C:\\Users\\james\\Documents\\DKI\\pyAFQ\\BIDS\\derivatives\\afq\\"

dwi_data_files = glob(infolder+"\\**\\ses-*\\*dwi*preprocessed.nii.gz")
bval_files = glob(infolder+"\\**\\ses-*\\*preprocessed.bval")
bvec_files = glob(infolder+"\\**\\ses-*\\*preprocessed.bvec")

dki_ak_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dki_ak.nii.gz")
dki_kfa_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dki_kfa.nii.gz")
dki_mk_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dki_mk.nii.gz")
dki_rk_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dki_rk.nii.gz")
dti_ad_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dti_ad.nii.gz")
dti_fa_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dti_fa.nii.gz")
dti_md_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dti_md.nii.gz")
dti_rd_scalar_files = glob(infolder+"\\**\\ses-*\\metrics\\*dti_rd.nii.gz")

output_dirs = []
for dwi_data_file in dwi_data_files:
    output_dir1 = dwi_data_file.replace("pydesigner", "afq")
    if len(output_dir1) == 113:
        output_dir2 = output_dir1.replace(output_dir1[72:127], "")
        output_dirs.append(output_dir2)
    if len(output_dir1) == 123:
        output_dir2 = output_dir1.replace(output_dir1[77:127], "")
        output_dirs.append(output_dir2)
    if len(output_dir1) == 127:
        output_dir2 = output_dir1.replace(output_dir1[79:127], "")
        output_dirs.append(output_dir2)

for x in range(0, len(dwi_data_files)):
    dki_ak_scalar_file = ImageFile(path = dki_ak_scalar_files[x])
    dki_kfa_scalar_file = ImageFile(path = dki_kfa_scalar_files[x])
    dki_mk_scalar_file = ImageFile(path = dki_mk_scalar_files[x])
    dki_rk_scalar_file = ImageFile(path = dki_rk_scalar_files[x])
    dti_ad_scalar_file = ImageFile(path = dti_ad_scalar_files[x])
    dti_fa_scalar_file = ImageFile(path = dti_fa_scalar_files[x])
    dti_md_scalar_file = ImageFile(path = dti_md_scalar_files[x])
    dti_rd_scalar_file = ImageFile(path = dti_rd_scalar_files[x])
    
    if not os.path.exists(output_dirs[x]):
        os.makedirs(output_dirs[x])
    
    myafq = ParticipantAFQ(
        dwi_data_file = dwi_data_files[x],
        bval_file = bval_files[x],
        bvec_file = bvec_files[x],
        output_dir = output_dirs[x],
        tracking_params=dict(
            odf_model="DTI"
            ),
        segmentation_params=dict(
            seg_algo="AFQ"
            ),
        scalars=[
            "dti_fa",
            "dti_md",
            "dki_fa",
            "dki_md",
            "dki_awf",
            "dki_mk",
            dki_ak_scalar_file,
            dki_kfa_scalar_file,
            dki_mk_scalar_file,
            dki_rk_scalar_file,
            dti_ad_scalar_file,
            dti_fa_scalar_file,
            dti_md_scalar_file,
            dti_rd_scalar_file
            ]
       )
    shutil.copy(r'C:\Users\james\Documents\DKI\pyAFQ\Template Images\tpl-MNI152NLin2009cAsym_res-01_desc-brain_mask.nii.gz',r'C:\Users\james\.cache\templateflow\tpl-MNI152NLin2009cAsym\tpl-MNI152NLin2009cAsym_res-01_desc-brain_mask.nii.gz')
    shutil.copy(r'C:\Users\james\Documents\DKI\pyAFQ\Template Images\tpl-MNI152NLin2009cAsym_res-01_T1w.nii.gz',r'C:\Users\james\.cache\templateflow\tpl-MNI152NLin2009cAsym\tpl-MNI152NLin2009cAsym_res-01_T1w.nii.gz')
    myafq.export_all()

The output I get is:

The visualizations accurately represent the lack of fibers found using v0.12.

It is not obvious to me what parameters may have been changed between versions to cause such a drastic change in the number of fibers tracked. I will note that I ran the stanford_hardi example and replicated the data on the pyAFQ website (AFQ API example).

Some assistance would be greatly appreciated.

[36000]

One major change in between those versions is that starting version 0.11 we changed tractography defaults from deterministic DTI to probabilistic CSD. In your case, if you were satisfied with version 0.8 outputs, you could try doing this:

    tracking_params=dict(
        odf_model="DKI",
        directions="det"
        ),

If the version 0.8 outputs also missed some bundles, another tactic is to increase number of seeds or to only seed in the ROI, like so:

    tracking_params={"n_seeds": 3,
                     "directions": "prob",
                     "odf_model": "CSD",
                     "seed_mask": RoiImage()},

this is from the optic radiations (OR) example, because the OR is harder to find: https://yeatmanlab.github.io/pyAFQ/auto_examples/plot_optic_radiations.html#sphx-glr-auto-examples-plot-optic-radiations-py

Of course, more seeds is more computationally expensive.

Edit: if re-running in the same folder, make sure to delete the .trk / .csv files before re-running so pyAFQ knows to recalc these.

[jdias001]

I took your advice, but I am running into the same problem. Even after changing odf_model=DTI (what I used with v0.8), directions="det", and n_seeds=3, I am still missing many of the tracts I was able to define in the earlier version of pyAFQ. Changing n_seeds will find more traces for the same tracts, but it will not find more tracts.

n_seeds=1

n_seeds=3

n_seeds=5

[36000]

Are you able to share the data or derivatives? If so I can debug on my end. If not, I think there are a few things to check. 1. Is the mapping correct? You can see this by seeing if ...template_xform.nii.gz mostly lines up with the dwi file (this is the MNI template transformed into subject space). 2. Are the ROIs in the right place? To check this either look at the individual bundle visualization in viz_bundles (even the bundles without streamlines found will show ROIs) or look at rois in the ROIS folder.

[jyeatman]

Or maybe the diffusion directions are flipped

…

Sent from my iPhone

________________________________ From: John Kruper ***@***.***> Sent: Friday, May 13, 2022 11:07:47 AM To: yeatmanlab/pyAFQ ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [yeatmanlab/pyAFQ] Substantially Fewer Steamlines and Tracts using v0.11/v0.12 compared to v0.8. (Issue #831) Are you able to share the data or derivatives? If so I can debug on my end. If not, I think there are a few things to check. 1. Is the mapping correct? You can see this by seeing if ...template_xform.nii.gz mostly lines up with the dwi file (this is the MNI template transformed into subject space). 2. Are the ROIs in the right place? To check this either look at the individual bundle visualization in viz_bundles (even the bundles without streamlines found will show ROIs) or look at rois in the ROIS folder. — Reply to this email directly, view it on GitHub<#831 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AANRRSWK7OQVTQZ6TWUBICDVJZV4HANCNFSM5VZHMPGQ>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[jdias001]

The template_xform.nii.gz overlays nicely with the preprocessed dwi data file.

The preprocessed dwi, bvec, and bval files are available here:
https://www.dropbox.com/s/uq4mse0mrukckf2/sub-1221_ses-Plasticity_dwi_preprocessed.zip?dl=0

[36000]

Is it true that you are using the same dwi files with ParticipantAFQ that GroupAFQ found?

Edit: One reason I ask this, is that I would not be surprised if GroupAFQ did not choose the preprocessed DWI file if it saw an un-preprocessed one that ends in _dwi.nii.gz

[jdias001]

Hmm... So this may be the problem. Running the same script on the raw, non-preprocessed data yields a lot more streamlines and tracts:

This now raises the question of what pydesigner is doing to the dwi.

[36000]

Very interesting!

[36000]

Let's leave this issue open until we understand this better

[arokem]

Hey @jdias001 : thanks for reporting this! My understanding is that you have stumbled on some kind of unfortunate interaction between the pydesigner preprocessing and DIPY's tractography (but I am not yet ruling out the possibility that pyAFQ is doing something wonky... I'd say we haven't nailed this down yet). Thanks for sharing those processed files. We'll take a look and see what's going on. If I had to guess, I'd say that something about pydesigner's processing is obscuring the directional information in the data, which prevents tractography from proceeding effectively. In the meanwhile, could you also share the raw data for the same subject/session that you shared the processed data for?

[jdias001]

Here is a link to a Dropbox folder containing the raw and preprocessed files:

https://www.dropbox.com/scl/fo/zos9ut5r7u404rdemw8jv/h?dl=0&rlkey=mte5it1t2izzeb27ip3ipbapy

[arokem]

Great. Thank you! I'll try to diagnose this and get back to you.

[arokem]

I don't know if this matters at all, but the processed data is RL flipped relative to the preprocessed data. The b-vectors are near-identical. The processed data also looks like it's been smoothed a bit. Here's the first b=0 volume (unprocessed followed by processed):

And the subtraction between the images (after flipping back so they're both in the same RL orientation):

[arokem]

OK. I think that I found the issue. It's probably the RL flip, together with the no flipping of the b-vecs. Here are the tensor ellipsoids from the corpus callosum with the raw data:

And with the pydesigner-processed data:

As @jyeatman correctly suggested above, the b-vectors in the processed data are oriented in the wrong directions relative to the data, making it impossible to track with these data.

(if you're curious how I made these images are created, it's based on this DIPY example: https://dipy.org/documentation/1.5.0/examples_built/reconst_dti/#example-reconst-dti)

[TheJaeger]

@arokem you are absolutely correct - PyDesigner flips the image into RAS orientation without flipping the b-vecs. We discovered this recently while attempting tractography with PyDesigner outputs. Needless to say, a fix is on the way.

[arokem]

Awesome! Glad to hear.

@jdias001 : I'll go ahead and close this issue, but please feel free to reopen if this does not address all of your concerns.