Custom tractography using pyAFQ and QSIRecon output

Hi everyone !

Summary of what happened:

I correctly pre-processed my diffusion data with QSIPrep and performed a reconstruction using QSIRecon following the following pipline:

mrtrix_multishell_msmt_pyafq_tractometry
https://qsirecon.readthedocs.io/en/latest/builtin_workflows.html#mrtrix-multishell-msmt-pyafq-tractometry

This pipeline works correctly and has enabled me to obtain tractograms of the main paths.

However, I have specific regions of interest that I’d like to exploit with certain rules (start, end, pass through). Using the examples available in the pyAFQ doc, I’ve tried to adapt my code to include my regions of interest while preserving the tractography generated by QSIRecon according to the piple chosen.

It seems to me that I should use the following file: sub-CTR01_ses-01_space-T1w_desc-preproc_desc-tracks_ifod2.tck

Here is my complete code:

import nibabel as nib
from AFQ.api.group import GroupAFQ
import AFQ.api.bundle_dict as abd
import os

# Load the existing tractogram
tractogram_path = '/scratch/lcorcos/EcriParkTracto/derivatives/qsirecon/sub-CTR01/ses-01/dwi/sub-CTR01_ses-01_space-T1w_desc-preproc_desc-tracks_ifod2.tck'

# Load the ROIs (in NIfTI format, subject space)
roi_base_path = "/home/lcorcos/EcriPark_Code/Custum_tracto/roi_sub-CRT01/"
rois = {
    "Left_PMd": roi_base_path + "Left_PMd.nii.gz",
    "Left_M1": roi_base_path + "Left_M1.nii.gz",
    "Left_S1": roi_base_path + "Left_S1.nii.gz",
    "Left_SMA": roi_base_path + "Left_SMA.nii.gz",
    "Left_IPL": roi_base_path + "Left_IPL.nii.gz",
    "Left_Putamen": roi_base_path + "Left_Putamen.nii.gz",
    "Left_Thalamus": roi_base_path + "Left_Thalamus.nii.gz",
    "Right_Anterior_Cerebellum": roi_base_path + "Right_Anterior_Cerebellum.nii.gz",
    "Right_Posterior_Cerebellum": roi_base_path + "Right_Posterior_Cerebellum.nii.gz",
}

# Define the bundle with constraints for filtering
bundle = abd.BundleDict({
    "Writing": {
        "start": [rois["Left_PMd"], rois["Left_M1"], rois["Left_S1"], rois["Left_SMA"]],

        "include": [rois["Left_IPL"], rois["Left_Putamen"], rois["Left_Thalamus"]],

        "end": [rois["Right_Anterior_Cerebellum"], rois["Right_Posterior_Cerebellum"]],

        "cross_midline": None,

        # Mahalanobis distance filtering parameters
        "mahal": {
            "clean_rounds": 10,
            "length_threshold": 4,
            "distance_threshold": 2
        }
    }
})

# Initialize the GroupAFQ object
my_afq = GroupAFQ(
    bids_path="/scratch/lcorcos/EcriParkTracto/",
    output_dir="/home/lcorcos/EcriPark_Code/Custum_tracto/outputV2",
    bundle_info=bundle,
    tracking_params={},
    import_tract={'scope': 'qsirecon', 'suffix': 'tracks_ifod2'}
)

# Export all results
my_afq.export_all()

# Filter the fibers that match the constraints
filtered_tractogram = my_afq.tractogram

# Save the filtered tractogram
nib.save(filtered_tractogram, os.path.join("/home/lcorcos/EcriPark_Code/Custum_tracto/output", "filtered_tractogram.trk"))

I’m using pyAFQ version 2.0:

The problem I’m having is importing and/or detecting the tractogram already produced. I’ve tried these two versions:

import_tract={'scope': 'qsirecon', 'suffix': 'tracks_ifod2'}

Which gives the following error: ValueError: No custom tractography found for subject CTR01 and session 01.

And

import_tract=tractogram_path

Which gives the following error: AttributeError: ‘list’ object has no attribute ‘get_fdata’

For testing purposes, I have isolated a subject, but always following the BIDS structure:

.
├── derivatives
│   ├── freesurfer
│   │   └── sub-CTR01
│   │       └── ses-01
│   │           ├── label
│   │           ├── mri
│   │           │   ├── orig
│   │           │   └── transforms
│   │           │       └── bak
│   │           ├── scripts
│   │           ├── stats
│   │           ├── surf
│   │           ├── tmp
│   │           ├── touch
│   │           └── trash
│   ├── qsiprep
│   │   └── sub-CTR01
│   │       ├── anat
│   │       ├── atlas
│   │       ├── figures
│   │       └── ses-01
│   │           ├── anat
│   │           └── dwi
│   └── qsirecon
│       └── sub-CTR01
│           ├── figures
│           └── ses-01
│               └── dwi
│                   └── sub-CTR01_ses-01_space-T1w_desc-preproc
│                       ├── bundles
│                       ├── clean_bundles
│                       ├── ROIs
│                       ├── viz_bundles
│                       └── viz_core_bundles
└── sub-CTR01
    └── ses-01
        ├── anat
        ├── dwi
        ├── fmap
        └── func

Have any of you already done something similar? If so, I’d love to hear from you!

Thank you very much,

Ludovic

Hi @MagicLudo,

Thanks for posting your question. Based on the second error you posted, I think that you might be hitting a limitation of the pyAFQ, which is that it doesn’t know how to deal with lists of ROI file-names as definitions of start/end/include for a tract. Instead, what is expected is a nibabel image object. So, I would suggest create an array that is a union of all the ROIs you would like to use as a start/end/include for the tract (e.g., [rois["Left_PMd"], rois["Left_M1"], rois["Left_S1"], rois["Left_SMA"]]) and then allocate that array into a nibabel Nifti1Image object, which you could then pass to the initializer of the BundleDict object.

Hope that helps!
Ariel