TinasheMTapera
commented
4 years ago Thanks for bringing this to our attention, I’ll get to investigating this issue ASAP
Closed bbiney1 closed 4 years ago
TinasheMTapera
commented
4 years ago Thanks for bringing this to our attention, I’ll get to investigating this issue ASAP
TinasheMTapera
commented
4 years ago Hey @bbiney1, you've brought up something interesting for us to discuss. I'll loop in @mattcieslak so he can also weigh in.
Basically, your strategy for incrementing runs isn't quite how fw-heudiconv works. In processing, the program grabs a list of sequence info objects and checks against your logic which to add to the naming template. If you have one session with 5 acquisitions, your run information will increment as expected. But if you have any more than one session, with your current logic you'll fill your run-1 naming template with whatever the program finds first, and then everything will else will be added to run-2, because both the sequence info list and the naming template list are static; they don't reset for each session.
In the past we've approached separating runs by simply nesting another if statement in the logic for a sequence, and further parsing out the sequence info object s. For example, BBL's effort task had two runs, but were easily separable by the protocol name:
| total_files_till_now | example_dcm_file | series_id | dcm_dir_name | dim1 | dim2 | dim3 | dim4 | TR | TE | protocol_name | is_motion_corrected | is_derived | series_description | sequence_name | image_type | series_uid |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 537 | 1.3.12.2.1107.5.2.43.167024.2018032214185313767709937.MR.dcm | 5c1a86f09011bd0013369952 | ep2d_effort1_236_6.nii.gz | 448 | 448 | 236 | -1 | 3 | 0.03 | ep2d_effort1_236 | False | False | ep2d_effort1_236 | _epfid2d1_64 | ('ORIGINAL', 'PRIMARY', 'M', 'ND', 'MOSAIC') | 1.3.12.2.1107.5.2.43.167024.2018032214135115699603908.0.0.0 |
| -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- |
| 771 | 1.3.12.2.1107.5.2.43.167024.2018032214343969862950630.MR.dcm | 5c1a86f09011bd001436a0e5 | ep2d_effort2_236_7.nii.gz | 448 | 448 | 234 | -1 | 3 | 0.03 | ep2d_effort2_236 | False | False | ep2d_effort2_236 | _epfid2d1_64 | ('ORIGINAL', 'PRIMARY', 'M', 'ND', 'MOSAIC') | 1.3.12.2.1107.5.2.43.167024.2018032214283117672136244.0.0.0 |
If you're able to parse your runs apart, you can do the same:
elif "restingbold" in protocol:
if "_mb6_1200" in protocol:
info[rest1].append(s.series_id)
else:
info[rest2].append(s.series_id)If not, we may have to look into some other strategy.
bbiney1
commented
4 years ago This doesn't apply to our data. Here's the tsv file for our project. This project has two runs of resting bold data.
From looking at this, it seems like the column for dcm_dir_name in the rows for the two rest bold runs have a number at the end. And it seems like the number for run 1 is smaller than the number for run 2.
But it does seem like the column "dcm_dir_name" might be able to distinguish runs 1 and 2. For run 1 this attribute is restingBOLD_mb6_1200_10.nii.gz For run 2 this attribute restingBOLD_mb6_1200_14.nii.gz
The number at the end is different. For run 1, there's a '10'. For run 2, there's a '14'
If we extract the numbers from the string in dcm_dir_name then try to compare the numbers, sorting them from smallest to largest, we can use that to distinguish which rows contain run 1 and run 2 (and hopefully build something that can scale for additional runs). Does this plan seem like it could work?
I'm going to write the code up for this right now, I apologize if this explanation didn't quite make sense.
TinasheMTapera
commented
4 years ago You're absolutely right, your protocols don't have much distinguishing different runs.
If we extract the numbers from the string in dcm_dir_name then try to compare the numbers, sorting them from smallest to largest, we can use that to distinguish which rows contain run 1 and run 2 (and hopefully build something that can scale for additional runs).
This strategy could certainly work as long as all of your dicoms are named appropriately for it. I know that we've had situations where the dicoms don't have consistent filenames, but hopefully yours are consistent. Let me know how it goes
bbiney1
commented
4 years ago We are still convinced that something weird is happening specific to running curate at the project level.
I made the following changes to the heuristic. We left the calls to get_both_series and the other helper functions created previously, those functions are defined in my original github issue post.
for s in seqinfo:
protocol = s.protocol_name.lower()
fileCount = s.total_files_till_now
# Baseline Anatomicals
if "t1w_mpr" in protocol:
get_series(t1w, s)
# Diffusion spectrum image
elif "dsi_1.8mm_257dir_b5000_mb3" in protocol:
get_series(dwi_257dir, s)
#Resting task scans (Need to change according to timepoint (date))
elif "restingbold_mb6_1200" in protocol:
dcm_dir_str = s.dcm_dir_name
ints_in_str = re.findall(r'\d+', dcm_dir_str)
identifier_int = int(ints_in_str[-1])
identifier_dict[s.series_id] = identifier_int
elif "b0map" in protocol and "M" in s.image_type:
get_multi_b0mag(fmap_run1_mag, fmap_run2_mag, s)
elif "b0map" in protocol and "P" in s.image_type:
get_both_series(fmap_run1_ph, fmap_run2_ph, s)
sorted_dict_lst = sorted(identifier_dict.items(), key=lambda kv: kv[1])
info[rest1].append(sorted_dict_lst[0][0])
info[rest2].append(sorted_dict_lst[1][0])
# https://stackoverflow.com/questions/613183/how-do-i-sort-a-dictionary-by-value
return infoI made the change I mentioned earlier about extracting the last digit from the dcm_dir_name column. For each extracted digit, I stored the digit as a value in a dictionary, identifier_dict. The keys in this dictionary are the corresponding series id. I sorted this dictionary from smallest to largest by its values, and selected the series id mapped onto the smaller extracted digit for run 1 and the series id mapped onto the larger digit for run 2.
So the behavior we got was... This heuristic successfully ran for 1 subject and 1 session. it produced run 1 and run 2 for not only rest scans, but also for the fieldmaps and so forth surprisingly (without me implementing any changes to those areas).
Upon switching to curating the entire project... It seems for every person that wasn't the 1 subject we first tested the new heuristic on, they only had data from run2 produced. For the 1 subject that we did test the new heuristic on previously, that 1 subject had run1 and run2 only for the resting scans. The run1 for the fieldmaps and other types of acquisitions disappeared. Even though running this heuristic on the whole project should have produced the same effects on every subject, rather than singling out the 1 subject we had previously tested the heuristic on. Also what this behavior indicates is that even with the changes I made to how we parse the resting scans, we still only get 'run2' for almost every subject.
So we're still a bit skeptical of how fw_heudiconv_curate is behaving on the project level.
TinasheMTapera
commented
4 years ago I don't think this is an issue with unexpected behaviour vis-à-vis whole project vs. single subject. As mentioned, fw-heudiconv collates a list of acquisitions and their sequence info that is static regardless of how big the initial query was. You can see this here; the default behaviour is that the initial query to the database collects the entire project's data, and the --subject --session flags are only used to filter that query. The rest of curation goes through a single pipeline, so there's no reason to expect the curation itself to execute anything differently based on these flags.
So we can narrow this down more, could you run a few tests for me at the command line?
fw-heudiconv-curate --project TNI --heuristic /Users/imaginglab/Downloads/heuristic_TNI_final.py --subject <SUBJECT_1>
fw-heudiconv-curate --project TNI --heuristic /Users/imaginglab/Downloads/heuristic_TNI_final.py --subject <SUBJECT_2>
fw-heudiconv-curate --project TNI --heuristic /Users/imaginglab/Downloads/heuristic_TNI_final.py --subject <SUBJECT_2> <SUBJECT_1>
In 1. and 2. you should notice correct behaviour, and for 3. you should only see one of them curated correctly.
bbiney1
commented
4 years ago I just want to make sure I'm getting on the same page in terms of understanding what's happening.
If the way we were trying to increment runs was wrong, then you're probably indicating that our helper functions for get_both_series() and get_multi_b0mag aren't going to work as intended.
I can see the potential for problems in those helper functions. It seems like right now, they're written to act as if the info dictionary is 'refreshing' or emptying itself for every session. For clarity, in infotodict() we aren't quite iterating over sessions. Rather, the for loop iterates over the seqinfo objects in the seqinfo parameter.
I want to assume that the seqinfo objects are ordered in seqinfo list such that you go over the seqinfo related to each type of scan acquisition for 1 subject/session, then you move on to another set of acquisitions from a different subject/session. Is this assumption about the ordering of the seqinfo list correct?
Regardless of my assumption... Once the for loop moves on to a seqinfo object for a new session after the 1st session it iterated over, problems occur. We will end up passing in the same parameters used in a previous session for key1 and key2, and only the second conditional will be true for everything except for data stored for the first session.
For instance, on the first session we should have no problem storing the correct session ids for rest1 and rest2 in the correct keys within the info dict. Then iterating over another session and checking that info[rest1] is not empty/has a length of more than 2 by the time we get to this next session, so we only ever end up actually populating info[rest2]).
the way I can only imagine that we have to overhaul the way we're storing acquisitions with more than one run. I've considered that one potential solution may be for me to basically reuse the code I wrote for extracting the last number from dcm_dir_name and sorting as needed. I am trying to think of simpler work arounds.
For clarity, what variables specifically represent the "sequence info list" and the "naming template list?" Does sequence info list refer to the seqinfo parameter? By "naming template list," are you referring to the values of the info dictionary for each key?
(I don't have flywheel access until Monday afternoon, but answering these questions would be really useful in order for me to help diagnose the problem)
TinasheMTapera
commented
4 years ago Hey @bbiney1, let me get right to your questions:
For clarity, what variables specifically represent the "sequence info list" and the "naming template list?" Does sequence info list refer to the seqinfo parameter? By "naming template list," are you referring to the values of the info dictionary for each key?
The "sequence info list" is a list of tuples generated from the dicom header data of all of the files in the user's query. In the code we define it as seq_infos. Unfortunately, this list is not ordered.
"Naming template list" refers to the formattable strings the user defines in the heuristic for their filenames, eg sub-{subject}_ses-{session}_task-object_viewing_run-01_day-01_bold.nii.gz.
It seems like right now, they're written to act as if the info dictionary is 'refreshing' or emptying itself for every session
Indeed, this is not the case. seq_infos is generated once at the time of querying for all of the data in your query, and the loop in your heuristic does not refresh it. As mentioned, most users have data in their dicom headers that they can use to parse out runs from one another; this is the first that we've come across that doesn't. I'm sorry this is proving to be difficult to work on.
We can consider he approach you describe, and look towards setting up a framework that will allow it. In the meantime, since you are running this from the command line, it might be worth your while to loop over each session as mentioned. You could do, for example:
for SES in sessions:
fw-heudiconv-curate --project TNI --heuristic /Users/imaginglab/Downloads/heuristic_TNI_final.py --session $SES
jaredpz
commented
4 years ago I have an idea for solving this problem, but being a python newby I'm not exactly sure code it. Let me try to explain my solution here.
I propose to make the infotodict function a two stage function where you first populate a dictionary that has every run of a given protocol_name assigned to a single dictionary key, but instead of just adding s.series_id to the dictionary, we add a tuple of [s.series_id, s.total_files_till_now]. Then in the second step, you split these dictionary elements up into separate dictionary keys for each corresponding run that you have and you can order them by total_file_till_now.
So in our current heuristic we're populating a dictionary that looks like this:
info = {
# t1
t1w: [],
# field map
fmap_run1_ph: [], fmap_run1_mag: [], fmap_run2_ph: [], fmap_run2_mag: [],
# rest
rest1: [], rest2: [],
# dwi
dwi_257dir: [],
}But I propose we first populate a dictionary that looks like this
info = {
# t1
t1w: [],
# field map
fmap_ph: [], fmap_mag: [],
# rest
rest: [],
# dwi
dwi_257dir: [],
}where fmap_ph: [], fmap_mag: [], and rest: [] all contain a list of tuples of the form [s.series_id, s.total_files_till_now]
Next, we would go through each element of this dictionary, and for any that are non-zero length lists we parse through the lists and split it up into different dictionary keys corresponding to run number and ordered based on total_file_till_now
@TinasheMTapera does this seem like a plan that would work? It feels like this approach (or something like it) would be a somewhat general heuristic structure that would be able to flexibly respond to variable scanning protocols with multiple runs of the same sequences.
For the record, when I originally wrote this heuristic I got the get_both_series() function from the example heuristic posted here: https://github.com/PennBBL/fw-heudiconv/blob/master/example_heuristics/multi-task_fmri.py
If that function is only going to work in particular contexts or only when running heudiconv on the session level it might be worth noting in a comment or something in the example heuristic file.
TinasheMTapera
commented
4 years ago Hey @bbiney1, we might have a fix for you here https://github.com/PennBBL/fw-heudiconv/blob/dc36e76d057ab9113afde3d42ecd11777a936955/fw_heudiconv/cli/curate.py#L84
In this PR 9c09a21, the seqinfo object is refreshed as it loops over each session. Would you like to try out v0.2.4_0.1.7?
@jaredpz we will definitely have to look at get_both_series too, it might not be a generalisable solution and hence shouldn't be on the example list.
bbiney1
commented
4 years ago Hi,
So if I'm interpreting this right, then this new version of curate would make it so that the seqinfos object is ordered by session?
From comparing versions of curate.py... It seems you've changed what you pass into get_seq_info() so that this function would no longer take the whole sessions variable as is. In the version of curate.py that we've been using, sessions is created by what's returned by the client.get_project_request() call. Then the entirety of sessions was passed into a function, get_seq_info().
But now this for-loop that you linked has been introduced, and it passes each individual session into get_seq_info() rather than the entirety of the sessions variable (with no further processing). And that I assume is going to actually make the seqinfo objects ordered/grouped by session.
Just want to be sure I'm explaining this correctly or mostly correctly.
TinasheMTapera
commented
4 years ago Not that seq_infos is ordered, but that it is now being generated only one session at a time, refreshing over each session. With this approach, you can expect the same result as when running fw-heudiconv-curate with a single argument in the --session flag. Can you give your original heuristic a try @bbiney1?
#!/usr/bin/env python
import os
import re
def create_key(template, outtype=('nii.gz',), annotation_classes=None):
if template is None or not template:
raise ValueError('Template must be a valid format string')
return template, outtype, annotation_classes
# Note: All files should be in nifti format before using this script
# Struct
t1w = create_key(
'sub-{subject}/{session}/anat/sub-{subject}_{session}_T1w')
# rest
rest1 = create_key(
'sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_run-1_bold')
rest2 = create_key(
'sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_run-2_bold')
# field maps
fmap_run1_ph = create_key(
'sub-{subject}/{session}/fmap/sub-{subject}_{session}_run-1_phasediff')
fmap_run1_mag = create_key(
'sub-{subject}/{session}/fmap/sub-{subject}_{session}_run-1_magnitude{item}')
fmap_run2_ph = create_key(
'sub-{subject}/{session}/fmap/sub-{subject}_{session}_run-2_phasediff')
fmap_run2_mag = create_key(
'sub-{subject}/{session}/fmap/sub-{subject}_{session}_run-2_magnitude{item}')
# dwi
dwi_257dir = create_key('sub-{subject}/{session}/dwi/sub-{subject}_{session}_dwi')
def infotodict(seqinfo):
"""Heuristic evaluator for determining which runs belong where
allowed template fields - follow python string module:
item: index within category
subject: participant id
seqitem: run number during scanning
subindex: sub index within group
session: session id
"""
info = {
# baseline
t1w: [],
# field map
fmap_run1_ph: [], fmap_run1_mag: [], fmap_run2_ph: [], fmap_run2_mag: [],
#tms session rest scans
rest1: [], rest2: [],
# dwi
dwi_257dir: [],
}
# sometimes patients struggle with a task the first time around (or something
# else goes wrong and often some tasks are repeated. This function accomodates
# the variable number of task runs
def get_both_series(key1, key2, s):
if len(info[key1]) == 0:
info[key1].append(s.series_id)
else:
info[key2].append(s.series_id)
def get_multi_b0mag(key1, key2, s):
if len(info[key1]) < 2:
info[key1].append(s.series_id)
else:
info[key2].append(s.series_id)
# this doesn't need to be a function but using it anyway for aesthetic symmetry
# with above function
def get_series(key, s):
info[key].append(s.series_id)
for s in seqinfo:
protocol = s.protocol_name.lower()
fileCount = s.total_files_till_now
# Baseline Anatomicals
if "t1w_mpr" in protocol:
get_series(t1w, s)
# Diffusion spectrum image
elif "dsi_1.8mm_257dir_b5000_mb3" in protocol:
get_series(dwi_257dir, s)
#Resting task scans (Need to change according to timepoint (date))
elif "restingbold_mb6_1200" in protocol:
get_both_series(rest1, rest2, s)
#elif "restingbold_mb6_1200" in protocol:
# info[rest].append(s.series_id)
elif "b0map" in protocol and "M" in s.image_type:
get_multi_b0mag(fmap_run1_mag, fmap_run2_mag, s)
elif "b0map" in protocol and "P" in s.image_type:
get_both_series(fmap_run1_ph, fmap_run2_ph, s)
return info
MetadataExtras = {
fmap_run1_ph: {
"EchoTime1": 0.00412,
"EchoTime2": 0.00658
},
fmap_run2_ph: {
"EchoTime1": 0.00412,
"EchoTime2": 0.00658
}
}
IntendedFor = {
fmap_run1_mag: ['{session}/func/sub-{subject}_{session}_task-rest_run-1_bold.nii.gz'],
fmap_run1_ph: ['{session}/func/sub-{subject}_{session}_task-rest_run-1_bold.nii.gz'],
fmap_run2_mag: ['{session}/func/sub-{subject}_{session}_task-rest_run-2_bold.nii.gz'],
fmap_run2_ph: ['{session}/func/sub-{subject}_{session}_task-rest_run-2_bold.nii.gz']
}
# Function to swap out the '.' in the subIDs as they are uploaded.
def ReplaceSubject(subj_label):
p = re.compile('\\.SC\\.1\\.')
return str(p.sub("", subj_label))Hi,
Yes it seems like this heuristic now works! We didn't run this again on our TNI project (we used bash scripting yesterday to curate each TNI subject individually), we ran it on a different project that contained functional data (but the heuristic it used was mostly copied and pasted from the heuristic for TNI, it uses the get_both_series() function). It seemed to work, things look right in BIDS view.
Thank you so much!
TinasheMTapera
commented
4 years ago Awesome stuff, glad it worked out!
Right now we're encountering behavior where, if we do the
fw_heudiconv_curateinstruction and pass in the additional --subject and --session flags for one subject, then that one particular subjects gets properly organized into the BIDS standard. We are able to see this reflected in the GUI and from usingfw_heudiconv_export.But if we run
fw_heudiconv_curatefor an entire project, if we don't use the --session and --subject flags... The GUI shows that upon trying to use BIDS view, no data from run 1 is included. We only get files from run 2 that are organized to BIDS standards. The files from run 1 just disappear.The print statements/log messages from
fw_heudiconv_curatedon't seem to be very informative. Here is a sample of our code.Please note that the print statements DO NOT CHANGE if I remove the ampersand.
Below is the code for a heuristic file. I don't think this should be helpful, but I"m including it just in case anyone is curious. Personally based on this code... I don't see how the conditionals for run2 could ever be activated without the conditionals for run1 occurring beforehand.