Qsiprep recon: have sift2 structural connectivity matrices been normalized by mu?

Software Support
1

Summary of what happened:

I ran qsiprep preprocessing (version 0.21.4), then reconstruction (version 0.19.1) via qsiprep --recon-only with --recon-spec mrtrix_singleshell_ss3t_ACT-hsvs. (Using singularity containers in both cases.) The reconstruction runs to completion (though no figures are generated; I understand that’s a problem in many HPC environments). I’m particularly interested in connectivity matrices generated by mrtrix3 (e.g., ...sift_radius2_countconnectome.csv found in the ...desc-exemplarbundles_dhollanderconnectome.zip archive file, or in, say, variable schaefer200_sift_radius2_count_connectivity found in the matlab file ..._desc-preproc_dhollanderconnectome.mat)

I would like to compare connectivity matrices across subjects. As I understand it, the tcksift2 program will generate a normalization factor “mu” (via -out_mu), specific to each subject, that can be used to adjust the connectivity matrices so that they can be compared across subjects.

[1] notes that this normalization factor is only valid when AFD intensity normalization has been performed in preprocessing. The only reference to intensity in the output of the qsiprep preprocessing phase alludes to the normalization of the B0 volumes across scanning sesssions, but I only have one such session per subject. So, I guess my questions are:

  • has the appropriate intensity normalization been applied in preprocessing?
  • in the case of the sift variants, has mu been applied to the connectivity matrices already?
  • if not, is mu stored somewhere for later use? (I haven’t found it as I search through the derivative files.)
  • finally, looking at equation 4 in [2], it looks like there might be enough info to regenerate mu from the files currently generated, but I’d rather not if possible

Any insight would be greatly appreciated — clearly, I’m not 100% sure of what I’m doing :frowning:

Many thanks,

  • Daniel

[1] Computing & comparing SIFT2-connectomes - MRtrix3 Community
[2] Smith RE, Tournier JD, Calamante F, Connelly A. SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography. Neuroimage. 2015 Oct 1;119:338-51. doi: 10.1016/j.neuroimage.2015.06.092. Epub 2015 Jul 8. PMID: 26163802.

 singularity run \
              --cleanenv \
              "${QSIPREP_SINGULARITY_IMAGE}" \
              "${inpt_bids_dir}" "${outp_dir}" participant \
              --verbose \
              --verbose \
              --stop-on-first-crash \
              --mem-mb "20000" \
              --recon-only \
              --recon-input "${outp_dir}/qsiprep" \
              --recon-spec "mrtrix_singleshell_ss3t_ACT-hsvs" \
              --fs-license-file "${FREESURFER_LICENSE}" \
              --freesurfer-input "${outp_func_dir}/sourcedata/freesurfer" \
              --participant-label "${subj}"

Version:

qsiprep 0.21.4 (for preprocessing)
qsiprep 0.19.1 (for reconstruction)

Environment (Docker, Singularity / Apptainer, custom installation):

Singularity

Data formatted according to a validatable standard? Please provide the output of the validator:

The data conforms to the BIDS standard

Relevant log outputs (up to 20 lines):

PASTE LOG OUTPUT HERE

Screenshots / relevant information:

2

Original poster here.

Perhaps I went into too much detail above. My basic question is, are the connectivity matrices computed via qsirecon meaningfully comparable across subjects? E.g., in the ultimate simple comparison, can I perform a two sample t-test between groups on, say, the connectivity/edge-value between node A and node B? Would there be any confounds I should include? (Intra-cranial volume estimates?)

If anyone can spare the time to confirm that this is possible, or what steps might be needed to make it possible, I would be very grateful!

Thanks!

  • Daniel
3

Hi @silophophe,

It depends on what the connectome is weighted by. SIFT2 to my understanding would not be comparable between subjects because different FRFs were used across subjects.

Best,
Steven

4

Understood (I think).

Reading the MRtrix3, there is a stage in their processing [1] where they

    1. normalize white matter intensity across subjects (via dwinormalize); and
    1. then generate an average response function (via responsemean) for the cohort that they use going forward in the tractography
      I’m guessing this is what you’re referring to?

I ran step 1 after qsiprep, but then ran into a dead end: is there an option to feed a predetermined, average response function into the qsirecon pipeline and have it bypass the individual subject estimations?

Thanks again for all your help! It’s an amazing project.

  • Daniel

[1] Fibre density and cross-section - Single-tissue CSD — MRtrix3 3.0 documentation

5

Hi @silophophe,

Correct, there’s no response averaging across subjects like you’d find in a fixed based analysis pipeline.

To my knowledge there’s no good way of incorporating your own files, and you might be best off right now doing all of the post processing separately.

Best,
Steven

6

Hi Steven,

Thank you for clarifying.

I can probably build my own pipeline based on MRtrix3 documentation. However, it would be really helpful if qsirecon output the actual commands being issued at each node of the pipeline. I’ve gone down to --verbatim x 3 (debug level) but I’m still not seeing actual shell commands. I can probably dig through some of the history stored in the mif files, but seeing the commands as they come up would be really helpful.

Is there any setting that can show me actual shell commands?

Thanks again.

Best,
Daniel

7

Hi @silophophe,

If you look around in the working directory you can look for command.txt files that show you the commands that were used. Since MRtrix is run in terminal, there should be command.txt files for each step (as opposed to processes that are completely Python-based). If you organize the node-specific folders by date added / date modified, you should get them roughly in order as they were applied.

I suppose you could also try to run this workflow independently and change it for your needs (qsiprep/qsiprep/workflows/recon/mrtrix.py at master · PennLINC/qsiprep · GitHub), but that would require some knowledge of Nipype. At the very least you can look at the functions called in this workflow for reference. That Gtihub file might also be good to cross-reference with what you see in the working directory. Sorry that there is no direct and intuitive answer for finding the commands…

Best,
Steven

8

Aha! Thank you so much! That’s awesome.

I’ve been submitting these (singularity) jobs to a compute node on a HPC. Either I haven’t specified a working directory, or I’ve set it up as local /tmp for performance considerations. In any event, the contents of the work dir are inaccessible that way. I’ll set one up on our /scratch directory so I have something to look at once the job finishes.

As for Nypipe, I know nothing about it, but having sorted through the its output so many times, I have a mental model of how it works — which may or may not reflect reality. But, as you suggested, just looking at the python file you linked, and using the command.txt files to see what branches were actually taken, I’m hopeful I can piece it all together.

No need to apologize; this is great. And thanks so much for your responsiveness! I think I’m set for now.

Best,
Daniel

9

Hi @silophophe,

Sounds good! Yes, please make sure you specify a working directory with the -w argument, and make sure it is mounted in singularity with the -B argument, e.g.:

singularity run \
              --cleanenv \
              -B $work:/work
              "${QSIPREP_SINGULARITY_IMAGE}" \
              "${inpt_bids_dir}" "${outp_dir}" participant \
              --verbose \
              --verbose \
              --stop-on-first-crash \
              --mem-mb "20000" \
              --recon-only \
              --recon-input "${outp_dir}/qsiprep" \
              --recon-spec "mrtrix_singleshell_ss3t_ACT-hsvs" \
              --fs-license-file "${FREESURFER_LICENSE}" \
              --freesurfer-input "${outp_func_dir}/sourcedata/freesurfer" \
              --participant-label "${subj}"
              -w /work

For what it its worth, I ran a fixel based analysis from QSIPrep derivatives separately without much issue. My code is here: GitHub - smeisler/Meisler_Reading_FBA: Code used in the Meisler and Gabrieli 2022 eLife paper on fixel based analyses of reading abilities.. It is meant to be run on HPCs (specifically with Slurm scheduler), maybe you’ll find it useful? In either case, if you find yourself stuck on anything feel free to reply here and I’ll do my best to help.

Best,
Steven

10

Steven,

This is fantastic, way more than I would have expected. Thank you!

This is a lot to digest; as you can probably tell, I’m not completely up to speed on diffusion weighted analysis. Let me sit with this a while and I may come back with additional questions once I’ve absorbed all this material.

I guess I do have one basic question: I’m beginning to grasp the significance of fixel-based analyses. As best I can understand, they are a way to meaningfully compare various dwi-related properties (e.g., fiber orientation/density/cross-section) at a voxel level. Aside from the use of a common tissue response for all subjects, there are subsequent steps in the fixel-based analysis that include, in particular, generating a group-level FOD template, warping each individual FOD to the template, and reorienting the FODs to account for any local change in direction due to the warp. I see how this would be necessary for a voxel-level comparison.

However, if I’m interested solely in “bulk” properties like connectivity between regions, is it sufficient(-ish) to use the group-derived response function for all subjects, as we discussed above? In your estimation, will implementing the full fixel-based-analysis pipeline make an appreciable difference to the tractography results?

Thank you for taking the time to help.

Best,
Daniel

11

There is a scaling factor you can calculate and apply to the sift weights so they can be compared across subjects. The mrtrix3_connectome bids app does this, but the qsiprep recon workflow stops short because there are currently no group-level workflows implemented.

12

Thank you @mattcieslak for weighing in.

Yeah! I had used the mrtrix3_connectome bids app before, whose third stage (“group”) informed the title of this thread. One of the reasons I switched to qsiprep is because I have old data with no reverse phase scans to compensate for susceptibility distortion. So instead, I used synb0-disco to generate an infinite-bandwith, distortionless B0 as the reverse phase component, and placed it in the fmap/ BIDS directory to trigger topup. At one point there was a fork of the mrtrix_connectome pipeline that relied on the developer version of MRtrix3 that allowed for feeding in the synb0 topup results directly into the lower-level eddy process. But that fork seems to have disappeared in the interim. So I’ve jumped to qsiprep, which has a lot of additional nice features.

Following the MRtrix3 suggested steps, after preprocessing with qsiprep, I could, in theory, compute the mean tissue response over my subjects and then, if qsirecon allowed for mean tissue response as an input, I could perform tractography with this mean response and the resulting connectivity matrices could be compared across subjects.

In the mrtrix3_connectome app, the adjustment by mu comes at the very end of the pipeline, after tractography. I gather, from what you are saying, that I can use the delta between individual and mean tissue response (computed before qsirecon) to come up with a correction for each subject (equivalent to mu?) that I could use to scale the qsirecon connectivity matrices so that they could be compared across subjects?

That could make my life much simpler, rather than reconstruct a full recon (hah!) pipeline.

Thanks so much for your input. I feel like I’m abusing the forum with so many (elementary?) questions, but I am grateful for your help.

Best,
Daniel

13

Hi @silophophe

There was a really nice paper that answers this better than I could here (Smith et al., 2022, Aperture Neuro): https://apertureneuro.org/article/81056-quantitative-streamlines-tractography-methods-and-inter-subject-normalisation. It talks about how/when connectomes are and aren’t comparable, and the various biases associated with different parameters. Hopefully it is digestible for someone not super familiar with DWI, but if not we can help clarify questions. A short answer to your question is that if you want to weigh your connectome by the estimated fiber densities, which serves to help correct false positives, you would need the SIFT scaling factor (which varies with the FRF) to be consistent across the cohort.

I don’t think using a subject-derived vs cohort averaged FRF for tractography would result in meaningully different streamlines, it is just the weighting that are applied to streamlines when making the connectome that really require that scaling factor / FRF to be consistent across participants.

Best,
Steven

14

Thanks Steve.

I’ll definitely go through the paper; I think I’m finally getting how all the pieces fit together, so maybe I’ll actually understand it.

Near term, I guess what I’m asking is: can I use the mean and and individual’s FRF to adjust that individual’s sift weights output by qsirecon so that the resulting connectivity matrices are comparable across subjects?

Or maybe put a different way: can I compute the equivalent of the mu correction factor generated by mrtrix_connectome by using the mean and individual FRF (that I can compute myself after qsiprep preprocessing)? As a naive guess, would mu = (average FRF) / (individual FRF)? I.e., if my FRF is below the average, I scale up the my sift weights accordingly, and vice-versa?

I’m sorry if I’m being hardheaded. Just trying to see if I can make use of the results output by qsirecon by adjusting after the fact, rather than redo the whole pipeline to incorporate the mean FRF at the outset.

Thanks!
Daniel

15

I looked through the mrtrix3_connectome code awhile back and it seemed like calculating mu wouldn’t be too hard. Then you can just rescale the weights. This should definitely go into the qsiprep documentation

16

Hi @mattcieslak,

That’s a good idea — I’ll check the mrtrix3_connectome code. The github readme for step 3 (group normalization) cites this preprint: OSF, and has a formula for mu. I should have cited that paper in my original post, too.

As long as the requisite inputs are stored in the qsiprep work directories, maybe this can work.

Thanks for your suggestion.

Best,
Daniel

17

Following up: I looked at source code for mrtrix3_connectome, and found the sections that compute the scale factor that normalizes the connectivity matrices so that they can be compared across subjects.

That scale factor is the product of four parts. Here’s the comment preceding the code that computes the scaling:`

# Third pass through session data in group analysis:
# - Scaling factors for connectome strengths:
#   - Multiply by SIFT proportionality coefficient mu
#   - Multiply by (mean median b=0) / (subject median b=0)
#   - Multiply by (subject RF size) / (mean RF size)
#     (needs to account for multi-shell data)
#   - Multiply by voxel volume

where RF is the tissue response.

I can get the B0- and tissue response-related values from mrtrix3 functions dwinormalize and responsemean. (Steps 5 and 6 in the MRtrix3 fixel-based analysis.) The voxel volumes are a given. That leaves factor mu from tcksift2. While mu doesn’t appear in the qsirecon’s dwi directory, looking at the mrtrix.py workflow in the qsiprep github repository, it does look like mu might get stored in the work directory:

`
…

        (tck_sift2, outputnode, [
            ('out_mu', 'mu'),
            ('out_weights', 'sift_weights')])

…
`

So I’ll run qsirecon with a work directory, verify that mu is present, and go from there.

Thank you all!
Daniel