Hello everyone,
I have already get the whole-brain tractogram with QSIRecon (mrtrix_multishell_msmt_ACT-hsvs). It has generated the track file (*model-ifod2_streamlines.tck,10 million) and the corresponding sift2 weight for each streamline (*model-sift2_streamlineweights.csv). Then I want to use pyAFQ to select some bundle of interest (BOI) by custom waypoint ROIs and probmaps.
According to my understanding, when calculating the metrics of a fiber bundle, the weighted average based on sift2 weight is required. So, the question is how can I calculate the scalar metrics (e.g. FA, MD) for each streamline and each bundle of my BOI with the streamline weight from sift2? Or, could it be possible to combine the sift2 weight and streamline together during the pyAFQ process of filtering the fiber bundles, and perform the filtering simultaneously?
That is a great question. I am not sure what a weighted average based on sift2 is “required”, but on the face of it, using sift2 weights seems like a great idea, and I would be quite interested to learn how/whether tract profiles that are weighted by sift2 weights differ from those that are weighted by the deviation of the streamline from the core of the tract (which is pyAFQ’s current default).
Fortunately, this should be possible to calculate with currently available software: The DIPY afq_profile function (see example here) takes an optional 'weights' key-word argument. I believe that if you provide the sift2 weights as an array that is ordered in the same order as your streamlines, the result would be a tract profile that is weighted proportionally by the sift2 weight assigned to each streamline.
I’d be quite interested to see a comparison if you conduct one.
Thank you very much for your kind reply. I agree that using the DIPY afq_profile function for weighted calculation of streamlines seems to make sense. But before that, we need to determine the index of the BOI streamlines in the original whole-brain tractogram, so that we can extract the sift2 weight in the right order. This is the problem that troubles me .
To achieve this, I have tried to write the original index into the metadata (data_per_streamline) of the whole-brain tractogram (both in .tck and .trk format) with DIPY. And then the tract file was provided to pyAFQ for bundle recognition, with the expectation that it can be filtered while carrying the index.
# Codes for writing the original index into the metadata
import numpy as np
from dipy.io.streamline import load_tractogram, save_tractogram
raw_file = "*desc-tracks_ifod2_100k.tck" # The same procedure applies to.trk files as well.
out_file = "*desc-tracks_ifod2_100k_withidx.tck"
sft = load_tractogramn(raw_file, reference="*space-T1w_FA.nii.gz"
n_streamlines = len(sft.streamlines)
orig_idx=np.arange(n_streamlines)
sft.data_per_streamline["orig_idx"] = orig_idx
save_tractogram(sft, out_file)
Unfortunately, it failed. I found that the .tck format does not seem to support metadata writing (this information is automatically erased during saving with DIPY save_tractogram function). Although the .trk format support metadata writing, it cannot be retained after pyAFQ bundle recognition…
Now I have another idea: For each streamline in the BOI, compare it one by one with the streamlines in the original whole-brain tractogram, thereby determining its index. But this seems to involve a relatively large amount of calculation (thousands of people * 10 million streamline). Do you think this plan is feasible? Or are there any other better alternatives?
By the way, although I had not thought about making a direct comparison between tract profiles that are weighted by sift2 weights and those under pyAFQ method before (I had only highly appreciated the bundle cleaning process of pyAFQ before), I now realize it is interesting. So after finding out the solution of the above problem, I’d like to do the comparison and share with you all
The recognize function, which is the work-horse of pyAFQ bundle recognition, has a return_idx key-word argument that is designed specifically for this kind of use-case, where indices in the original tractogram need to be propagated. This does mean that you will have to set up everything and call it directly rather than depending on pyAFQ’s data management and workflow API. Does that provide a good enough solution?
Correction: I think that I might be wrong, and you can actually pass return_idx=True to the GroupAFQ class and also set that up as part of the config file, if you are using the CLI. I believe that this would save the original indices in the individual tract trk files, just like what you were trying to do. So no need to write all the low-level code
Yes, just as you noticed, I also found that the clean function also has a return_idx parameter.
return_idx: bool
Whether to return indices in the original streamlines. Default: False.
In general, the cleaning process needs to be completed through multiple iterations (defaults to 5 iterations in pyAFQ). So the indices recorded in the log file seem to be relative to the most recent tractogram iteration, as shown in the following screenshots.
It seems that some modifications to the code of pyAFQ are necessary (although for me, a novice in coding, this might take some time). Nevertheless, this is indeed a promising discovery.
.
