automaticanalysis / automaticanalysis

Automatic Analysis (aa)
http://automaticanalysis.github.io
MIT License
72 stars 35 forks source link

QUESTION - How to deal with more than one anatomical scan in BIDS format? #331

Closed TamerGezici closed 2 years ago

TamerGezici commented 2 years ago

Dear all,

I am still new to BIDS and AA so if the answer to my question is out there somewhere, my apologies.

Let's say during a single testing session the participant undergoes fMRI and after a couple of runs they ask to be removed from the scanner (for a short toilet break or some discomfort they had to due to the scanner). When they are ready again, we put them back into the scanner and take a new anatomical image and other BOLD scans. Assuming my data set is BIDS compliant, how should I map such thing to BIDS so that AA will understand which T1 image to use?

And let us assume the paradigm / experiment that was being tested prior to the break is different than the one after the break. How can I map things so that AA will know these T1's will be used for seperate kinds of BOLDs?

Thanks.

jones-michael-s commented 2 years ago

If I am understanding the issue correctly, I believe BIDS "session" selection may provide a solution.

As an example, consider the ds000114 dataset (available on OpenNeuro and also used in a couple of the aa example scripts).

This dataset includes two sessions (these happen to be test and retest sessions done on different days, but they could have been on the same day before and after a break as you describe in your situation).

Here is how the data is structured in the BIDS data for a typical subject:

image

Data from the test and retest sessions are organized into the ses-test and ses-retest folders. Each of these has its own anatomical and functional data (as well as dwi, which is not relevant to your example).

In aa, you can select which session data to analyze by passing the session name to aas_processBIDS:

aap = aas_processBIDS(aap,{'sess-test'});

or

aap = aas_processBIDS(aap,{'sess-retest'});

The one difference is you mention using a different experimental paradigm before/after the break. I have never tried this in a multi-session BIDS dataset, but I can't think of an a priori reason why it would be a problem.

Hope this helps, Mike

tiborauer commented 2 years ago

Just a follow-up: According to the BIDS specification, "session" usually refers to visits, and "if a subject has to leave the scanner room and then be re-positioned on the scanner bed, the set of MRI acquisitions will still be considered as a session". Therefore, your use case seems to be a single-session one to me.

You can also have multiple structural scans within a single session, and you consider and model them as "run"s in BIDS: sub-01_run-1_T1w.nii.gz sub-01_run-2_T1w.nii.gz (only nonnegative integers are allowed as run labels)

If you have multiple structural acquisitions (which you do), you have to tell aa which to use. You do this by setting aap.options.autoidentifystructural_chooselast = 1; or aap.options.autoidentifystructural_choosefirst = 1; in a mutually exclusive way.

TamerGezici commented 2 years ago

Thank you for your response. I have a different question now.

My PI asked me to put our data in BIDS format (since it would be necessary for publications) but not to analyze the data using the process bids way (for more flexibility).

Do you have any example scripts to analyze data in BIDS format but without following the regular bids process method? I have sample scripts of him which he used in CBU which directly takes DICOM outputs. But it looks like CBU's way of processing things is different than what is supposed to be done for a BIDS dataset.

For example, just for testing purposes, I tried to define a subject like this (the code I am sharing here is simplified)

aap.acq_details.root = 'E:\fMRI\fr-pr-3-AA\HCFprep\AAresults';
aap.directory_conventions.rawdatadir = 'E:\fMRI\fr-pr-3-AA\fmri_raw_bids\bids';
aap.directory_conventions.analysisid = 'AA';
aap.directory_conventions.subject_directory_format = 1;

aap = aas_addsubject(aap,'sub-01','structural',{'E:\fMRI\fr-pr-3-AA\fmri_raw_bids\bids\sub-01\ses-1\anat\sub-01_ses-1_run-1_T1w.nii'});%28
aap = aas_addsubject(aap,'sub-01','functional',{'E:\fMRI\fr-pr-3-AA\fmri_raw_bids\bids\sub-01\ses-1\func\sub-01_ses-1_task-discrB_run-01_bold.nii'});%28

Error: Expected one output from a curly brace or dot indexing expression, but there were 2 results.

I inspected the aas_addsubject.m file to understand what all I should feed into this function, but after several attempts I couldn't get anywhere.

The default parameter file I edited based on looking at the CBU parameters:

<?xml version="1.0" encoding="utf-8"?>
<aap xmlns:xi="http://www.w3.org/2001/XInclude">
   <xi:include href="E:\fMRI\automaticanalysis\aa_parametersets\aap_parameters_defaults.xml"
               parse="xml"/>
   <local>
      <directory_conventions>
         <fieldmapsdirname desc='Subdirectory in subject dir for fieldmaps' ui='dir'>fmap</fieldmapsdirname>
         <structdirname desc='Subdirectory in subject dir for MRI structurals' ui='dir'>anat</structdirname>
         <specialseriesdirname desc='Subdirectory in subject dir for special series' ui='text'>func</specialseriesdirname>
         <rawdatadir desc="Root on local machine for processed data" ui="dir">E:\fMRI\fr-pr-3-AA\fmri_raw_bids\bids</rawdatadir>
         <eventsdirname  desc='Subdirectory in study, subject and/or session for event files' ui='text'>events</eventsdirname>
         <subjectoutputformat desc='sprintf formatting string to get subject directory from number' ui='text'>sub-%03d*</subjectoutputformat>
         <subject_directory_format desc='Format of subject directory' ui='optionlist' options='from subject_directory_names|from data|S#|manual'>2</subject_directory_format>
         <protocol_structural ui='text'>T1</protocol_structural>
         <toolbox desc="Toolbox with implemented interface in extrafunctions/toolboxes"
                  ui="custom">
            <name desc="Name corresponding to the name of the interface without the &#34;Class&#34; suffix"
                  ui="text">spm</name>
            <dir ui="dir_list">C:\Users\tamer\OneDrive\Documents\MATLAB\spm12</dir>
         </toolbox>
      </directory_conventions>
              <options>
            <wheretoprocess desc='where to do processing' options='localsingle|localparallel|aws|qsub' ui='optionlist'>qsub</wheretoprocess>
        </options>
      <acq_details>
         <root desc="Root on local machine for processed data" ui="dir">E:\fMRI\fr-pr-3-AA\HCFprep\AAresults</root>
      </acq_details>
   </local>
</aap>
tiborauer commented 2 years ago

I am not sure what you mean by "not to analyze the data using the process bids way". BIDS is about how to store data rather than about how to analyse it. With aa, you can be as flexible with your analysis as you want.

If your data is in BIDS format, then you can use the aas_processBIDS function as you can see in the tutorial or in the demos ds114 and ds114motor.

In case you wonder about the function aas_processBIDS, it does not do any actual processing (i.e. data manipulation) or analysis. It simply reads in the BIDS dataset into aa. Then, it is up to you how you set up your processing and analysis workflow.

TamerGezici commented 2 years ago

Ah I see. Yeah, I checked the aap structure after processing BIDS. I guess I can work around with this as well. That makes sense.

I asked that because we don't want to bother ourselves for now during our analysis by putting our events in BIDS format, using .json's and such. So I wanted clarification to see if that would be possible. For example, feeding event files on our own way to AA to do 1st and 2nd level analysis. I can see that this is pretty possible, by just working around with some scripting it should be OK.

I have another set of more questions, in terms of what's going on, what is the main difference between the preprocessing steps of SPM_CH30.xml and fmri.xml?

I can see that it takes dicom input, that's there, but I am more interested in the difference between coreg_extended_1 and 2, versus coreg_extended (in SPM_CH30).

For example, what does reorienttomiddle epi do? That is not there in the fmri.xml file. Also, slicetiming is explicitly defined in fmri.xml, but it is not in SPM_CH30, so which module is responsible for slice timing in the tutorial xml?

I know I am asking a lot of questions, but I'm still new to these as well, and I want to understand what I'm working with.

Thanks a lot for your help so far.

tiborauer commented 2 years ago

The main difference between the two pipelines is that the fmri processes fieldmaps and, consequently, uses realignunwarp instead realign. The one- and two-stage approaches for coreg_extended do not make a difference here; and the two-stage approach is more useful, if coregistration is required in more than one stage, e.g. if you have multiple modalities (e.g. fMRI, DWI) or more complex workflow (e.g. modelling, decoding).

TamerGezici commented 2 years ago

I see, thanks. And what does reorienttomiddle do?

tiborauer commented 2 years ago

It corrects for the difference in indexing for FSL and SPM. For FSL (and BIDS nii.gz), indexing starts at the corner, while for SPM, indexing starts in the middle of the brain (AC for MNI). The corner-based indexing confuses SPM registration.

TamerGezici commented 2 years ago

While running a preprocessing pipeline myself, I ran into Specified TPM $SPMDIR/tpm/TPM.nii not found. error.

But the file is there, in the correct directory, so there should not be an issue.

cfg.tpm = fullfile(spm('dir'), 'tpm', 'TPM.nii');

I added this to line 91 in aamod_segment8.m as a hacky solution.

cfg.tpm was assigned this:

image

I'm guessing it's an OS issue.

tiborauer commented 2 years ago

I have to take a look. It seems the $SPMDIR environmental variable is not evaluated.

You can overwrite the TPM specification in a non-hacky way in the user script: aap.tasksettings.aamod_segment8.tpm = fullfile(spm('dir'), 'tpm', 'TPM.nii');

tiborauer commented 2 years ago

Yes. As suggested the $SPMDIR environmental variable is not evaluated, because we have not implemented a straightforward way for expanding environmental variables in Windows (see #289).

TamerGezici commented 2 years ago

I see. Thanks for the quick fix.

I found a bug, and implemented a quick fix myself to the source code. The explanation may be a bit complicated, so I will share the fix first. The code to be fixed is this in aamod_epifromnifti.m, line 19:

I had to replace "1" in the original source with "sess".

        if ~iscell(series) ... 
                || (~isstruct(series{sess}) ... % hdr+fname
                && ~ischar(series{sess}) ... % fname
                && ~iscell(series{sess})) % fname (4D)
            aas_log(aap,true,['ERROR: Was expecting list of struct(s) of fname+hdr or fname in cell array\n' help('aas_addsubject')]);            
        end

Reproduction:

While the script was running for subject 3, I ran into the error message above.

I have plenty of different tasks in my experiment. A list of them is below.

image

Subject 1 and 2 were discrB subjects, subject 3 was a discrA subject. As far as I understand, the series variable holds a seperate index for each task type. My participants so far until subject 3, were on the task type discrB, which was assigned index 1 in sessions, so the code worked fine (although by chance). But the moment the task index to be processed was other than 1, it crashed.

So the code should be indexing the session number here, not 1, because then it won't be covering the corresponding task for that participant, right?

For example, look at subject 4's series. This participant also had done discrA:

image

If series is indexed at index 1, then the header and file name will never be found, because the code is looking for it in the wrong index.

At least that's what I understood, and implementing this fix made it work just fine.

tiborauer commented 2 years ago

Thank you for being so patient and thorough! :)

The use case you describe is a sparse one, i.e. subjects have different sets of tasks. I wonder whether it is really the best representation of your data. What is the difference between discrA, discrB, and discrC? If the acquisition parameters are the same, then you can use the same JSON to describe them. In that case, all subjects would have the same sets of tasks, and the code would work as was.

tiborauer commented 2 years ago

Your fix would still make sense, of course. Would you care sending a PR?

TamerGezici commented 2 years ago

Sure yes I can send a PR.

Discr A, B, C are similiar experiments which run on different types of blocks (the letters here signify different blocks). A block is a rule switch experiment (2 rules), B also has 2 rules but the task flow is hierarchically organized in a different way, and block C involves 4 rules and is also hierarchically organized. Easydiff is the same experiment but also hierachically organized in a different way than all others, and easydiff also has its own versions for each block (easydiffA,B,C and so on)

So they are all the same experiment but with different conditions. The experiments' structure are mainly similar, but there are subtle differences between them.

If what is meant by acquisition parameters of these experiments is the BOLD sequences, yes. The only difference is the behavioural parts of things. ie, the scripts executed are different.

I'll need to look into how these tasks can be made into a single .json I haven't seen any examples of something like that in openneuro yet. This was what made sense for me and I understood. In the end, all these are different tasks (at least in my mind) because the experiments do differ from each other in several ways.

TamerGezici commented 2 years ago

I ran my analysis for 20 participants and it ran fine. Thanks a lot. I have lots of disk space, but some members in my lab do not. We will be mostly using swar and war files for our analysis. Is there a way to delete the module outputs as the script runs in upcoming steps? Because I can see that the files will quickly take lots of space as the analysis is going on, and if I can't delete them during the process, the analysis will have to end at some point.

I remember seeing some setting somewhere, but I dont remember what it was.

tiborauer commented 2 years ago

First, I suggest opening a new issue for every independent issue :) you come across. This would make the individual threads more manageable.

There is an option in aa (aap.options.garbagecollection), which - if enabled - clears the duplicates. That already saves a lot of space.

If you are absolutely sure that you do not need intermediate derived data (e.g. for QA), then you can delete them. However, if you delete the diagnostic images (jpgs), then the report will not be able to display them, because they are linked rather than embedded into the HTMLs. Also deleting of the doneflags (done_* files) would trigger the re-execution of the given and all subsequent modules in the workflow in case of rerunning. It is not an issue if you are happy with the results and have no plan to rerun or extend the workflow.

TamerGezici commented 2 years ago

Thank you for your answer. Sure I will create a new issue for each independent one. I thought it would clutter less this way, but the other way would also be better for future reference.

TamerGezici commented 2 years ago

Dear all,

I am re-opening this for clarification.

I am using addsubject now to add subjects to AA (not processbids).

I am not sure how utilizing aap.options.autoidentifystructural_chooselast = 1; or aap.options.autoidentifystructural_choosefirst = 1; can help if a participant undergoes some amount of fMRI runs, leaves the scanner, enters it again and does more tasks. I need to use both anatomical scans.

Is there a way to specify in code to say (in pseudo code): Use these series with this anatomical image, and these series with this anatomical image.

tiborauer commented 2 years ago

First, you have to decide whether and how you want to use both anatomical scans.

TamerGezici commented 2 years ago

Thanks @tiborauer

The solution I had on top of my mind was also to use your last suggestion. But I don't know how straightforward that can be. I think I'd need to merge the folders together without touching the "done" flags in the smoothened files and continue with my GLM, right?

I am surprised that the realignment algorithm can handle even after the participant is re-positioned. I did some research but couldn't find anything on the capabilities of the realignment algorithm and how it can handle slight misalignments. How much would you consider "slight"? I feel like patient heads are usually very differently positioned than their previous anatomical scans (we usually let them completely leave the scanner for a bathroom break if they complain, or if they complain that their head hurts due to the cushions)

It may be worth running an experiment on my own and to see how much the results change when one of the anatomical scan is used.

tiborauer commented 2 years ago

I have good experience with realignment between repositionings. Of course, we also pay attention that the FOV of the acquisition should correspond to the new head position. The larger the misalignments, the less successful the realignment is. There are also special sequences, such as the AutoAlign from Siemens, which perform very well in repositioning. The largest misalignment between re-positioning (without AutoAlign) I encountered was <1cm. Ultimately, however, "you should always look at your data" (a key lesson learned during the FSL course). Try realignment across all sessions regardless of repositioning. If the output looks good, then carry on. If not try other approaches.

Registering the two positionings as separate subjects can be a straightforward solution, however, it can complicate the second-level design.

Another, not-yet-implemented, approach is to keep all data in one subject but realigning sessions separately and then coregistering them via their meanepis.

TamerGezici commented 2 years ago

Of course, we also pay attention that the FOV of the acquisition should correspond to the new head position

To check whether the previous FOV covers the new position of the head, would taking a quick localizer and inspecting how well it covers the head be a good solution? And if you see huge misalignments, taking a new anatomical scan? And is there any lee-way for the FOV? Could the realignment algorithm figure out to re-adjust the FOV on its own if the misalignment is small enough?

I think AutoAlign is not an option for us because we use a Tim Trio scanner. Do you know if Tim Trio supports it? I googled it and found this in siemens website. Are they the same? Our software version is B17, so it does match. Maybe I should ask this to the engineer in our facility, though :)

tiborauer commented 2 years ago

To check whether the previous FOV covers the new position of the head, would taking a quick localizer and inspecting how well it covers the head be a good solution? And if you see huge misalignments, taking a new anatomical scan?

Yes

And is there any lee-way for the FOV? Could the realignment algorithm figure out to re-adjust the FOV on its own if the misalignment is small enough?

No. Realignment is post-acquisition, therefore, it cannot change the FOV. Actually, it reduces the effective FOV, because it takes the overlap between the two (or more) FOVs. What you described is what AutoAlign does, which is pre-acquisition.

I have used AutoAlign with TIM Trio, when I was in Göttingen and also at the MRC CBU. However, check with your Siemens engineer.

TamerGezici commented 2 years ago

Thank you for your answers. This is very enlightening and eye-opening information.

Try realignment across all sessions regardless of repositioning.

I didn't get this part. Can you elaborate, please? My preprocessing protocol already does realignment (I guess) for all sessions. Is there something I need to do apart from adding the necessary module to my .xml file?

I have been using aamod_realign, shall I instead opt in for aamod_realignunwarp to utilize these things?

EDIT: I see. The latter requires a fieldmap. Well I will be using fieldmaps later on. I guess I cant use this for now.

tiborauer commented 2 years ago

I didn't get this part. Can you elaborate, please? My preprocessing protocol already does realignment (I guess) for all sessions. Is there something I need to do apart from adding the necessary module to my .xml file?

No. That is all. You need to check the output including the images and the moco plots.

I have been using aamod_realign, shall I instead opt in for aamod_realignunwarp to utilize these things?

EDIT: I see. The latter requires a fieldmap. Well I will be using fieldmaps later on. I guess I cant use this for now.

From your POV, aamod_realign and aamod_realignunwarp does not differ. And, indeed, if you have fieldmaps, you should probably use them (in aamod_realignunwarp).

TamerGezici commented 2 years ago

Thank you @tiborauer for your detailed and helpful answers. If I need to ask anything related to this topic again I will re-open the issue. For now, this is all I needed to know so I can close this issue.