Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

dwi preprocessing using metadata #15

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open

dwi preprocessing using metadata #15

Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
d97f843
shorter task label name
yeunkim Dec 2, 2016
ce1f242
dwi processing (multiple directions, runs)
yeunkim Dec 2, 2016
92cde71
run label change
yeunkim Dec 2, 2016
d7823d6
output directory name change
yeunkim Dec 5, 2016
6037f36
uncomment diff function
yeunkim Dec 5, 2016
31b460b
bug fixed
yeunkim Jan 26, 2017
6c256db
changed Dockerfile (npm update bug fixed)
Jan 27, 2017
86926c2
fixed repetitive for loops
yeunkim Mar 22, 2017
15ed837
Remove .idea from repo
Apr 11, 2017
0b41059
test drive
yeunkim Apr 16, 2017
6aa43ed
test drive
yeunkim Apr 16, 2017
db56702
Merge branch 'dwi' of https://github.com/yeunkim/HCPPipelines into dwi
yeunkim Apr 18, 2017
0eae829
uses metadata for processing
yeunkim Apr 20, 2017
00c3891
Merge branch 'master' of https://github.com/BIDS-Apps/HCPPipelines
yeunkim Apr 25, 2017
b79f6a8
Merge branch 'dwi_patch'
yeunkim Apr 25, 2017
2317eef
added Chris's suggestions and added more comments
May 3, 2017
4f6f47c
misspelling in comment
May 3, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions Dockerfile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,10 @@ ENV PATH /opt/freesurfer/bin:/opt/freesurfer/fsfast/bin:/opt/freesurfer/tktools:
# Install FSL 5.0.9
RUN apt-get update && \
apt-get install -y --no-install-recommends curl && \
curl -sSL http://neuro.debian.net/lists/trusty.us-ca.full >> /etc/apt/sources.list.d/neurodebian.sources.list && \
curl -sSL http://neuro.debian.net/lists/xenial.us-ca.full >> /etc/apt/sources.list.d/neurodebian.sources.list && \
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did you need to change the base Ubuntu distribution?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm. I actually wasn't aware of this change. There wasn't any need - sorry about that.

apt-key adv --recv-keys --keyserver hkp://pgp.mit.edu:80 0xA5D32F012649A5A9 && \
apt-get update && \
apt-get install -y fsl-core=5.0.9-3~nd14.04+1 && \
apt-get install -y fsl-core=5.0.9-1~nd+1+nd16.04+1 && \
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are those changes still necessary applicable to this PR? (applies to connectome-workbench as well)

rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*

# Configure environment
Expand All @@ -72,7 +72,7 @@ ENV FSLOUTPUTTYPE=NIFTI_GZ
RUN echo "cHJpbnRmICJrcnp5c3p0b2YuZ29yZ29sZXdza2lAZ21haWwuY29tXG41MTcyXG4gKkN2dW12RVYzelRmZ1xuRlM1Si8yYzFhZ2c0RVxuIiA+IC9vcHQvZnJlZXN1cmZlci9saWNlbnNlLnR4dAo=" | base64 -d | sh

# Install Connectome Workbench
RUN apt-get update && apt-get -y install connectome-workbench=1.2.3-1~nd14.04+1
RUN apt-get update && apt-get -y install connectome-workbench=1.2.3-1~nd16.04+1

ENV CARET7DIR=/usr/bin

Expand Down
111 changes: 81 additions & 30 deletions run.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,7 @@ def run_diffusion_processsing(**args):
'--echospacing="{echospacing}" '+ \
'--PEdir={PEdir} ' + \
'--gdcoeffs="NONE" ' + \
'--dwiname="{dwiname}" ' + \
'--printcom=""'
cmd = cmd.format(**args)
run(cmd, cwd=args["path"], env={"OMP_NUM_THREADS": str(args["n_cpus"])})
Expand Down Expand Up @@ -193,7 +194,7 @@ def run_diffusion_processsing(**args):
'fMRISurface', 'DiffusionPreprocessing'],
default=['PreFreeSurfer', 'FreeSurfer', 'PostFreeSurfer',
'fMRIVolume', 'fMRISurface',
'DiffusionPreprocessing'])
'DiffusionPreprocessing', 'TaskfMRIAnalysis'])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Leftover from another PR? This one should be just about DWI right?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, yes.

parser.add_argument('--license_key', help='FreeSurfer license key - letters and numbers after "*" in the email you received after registration. To register (for free) visit https://surfer.nmr.mgh.harvard.edu/registration.html',
required=True)
parser.add_argument('-v', '--version', action='version',
Expand Down Expand Up @@ -343,14 +344,15 @@ def run_diffusion_processsing(**args):
type='bold',
extensions=["nii.gz", "nii"])]
for fmritcs in bolds:
fmriname = "_".join(fmritcs.split("sub-")[-1].split("_")[1:]).split(".")[0]
fmriname = "_".join(fmritcs.split("sub-")[-1].split("_")[1:-1]).split(".")[0]
assert fmriname

fmriscout = fmritcs.replace("_bold", "_sbref")
if not os.path.exists(fmriscout):
fmriscout = "NONE"

fieldmap_set = layout.get_fieldmap(fmritcs)
print(fieldmap_set)
if fieldmap_set and fieldmap_set["type"] == "epi":
SEPhaseNeg = None
SEPhasePos = None
Expand Down Expand Up @@ -409,31 +411,80 @@ def run_diffusion_processsing(**args):
if stage in args.stages:
stage_func()

dwis = layout.get(subject=subject_label, type='dwi',
extensions=["nii.gz", "nii"])

# print(dwis)
# acqs = set(layout.get(target='acquisition', return_type='id',
# subject=subject_label, type='dwi',
# extensions=["nii.gz", "nii"]))
# print(acqs)
# posData = []
# negData = []
# for acq in acqs:
# pos = "EMPTY"
# neg = "EMPTY"
# dwis = layout.get(subject=subject_label,
# type='dwi', acquisition=acq,
# extensions=["nii.gz", "nii"])
# assert len(dwis) <= 2
# for dwi in dwis:
# dwi = dwi.filename
# if "-" in layout.get_metadata(dwi)["PhaseEncodingDirection"]:
# neg = dwi
# else:
# pos = dwi
# posData.append(pos)
# negData.append(neg)
#
# print(negData)
# print(posData)
posData = []
negData = []
PEdir = "None"
dwiname = "Diffusion"
dirnums = []
onerun = False

numruns = set(layout.get(target='run', return_type='id',
subject=subject_label, type='dwi',
extensions=["nii.gz", "nii"]))
# accounts for multiple runs, number of directions, and phase encoding directions
if not numruns:
onerun= True
numruns = {'run-01'}
if numruns:
for session in numruns:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for run in numruns ?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

run is a function - I believe. Would you like me to rename the variable session?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this would make more sense since you are iterating over runs not sessions.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I could change it to numrun ?

if not onerun:
bvals = [f.filename for f in layout.get(subject=subject_label,
type='dwi', run=session,
extensions=["bval"])]
else:
bvals = [f.filename for f in layout.get(subject=subject_label,
type='dwi', extensions=["bval"])]
dwi_dict = {'bvalFile':[], 'bval':[], 'dwiFile':[], 'direction':[]}
for bvalfile in bvals:
with open(bvalfile) as f: # get number of directions
bvalues = [bvalue for line in f for bvalue in line.split()]
dwi_dict['bvalFile'].append(bvalfile)
dwi_dict['bval'].append(len(bvalues) - 1)
dwiFile = glob(os.path.join(os.path.dirname(bvalfile),'{0}.nii*'.format(os.path.basename(bvalfile).split('.')[0]))) # ensures bval file has same name as dwi file
assert len(dwiFile) == 1
dwi_dict['dwiFile'].append(dwiFile[0])
dwi_dict['direction'].append(layout.get_metadata(dwiFile[0])["PhaseEncodingDirection"][0])

# check if length of lists in dictionary are the same
n = len(dwi_dict['bvalFile'])
assert all(len(dwi_dict[k]) for k,v in dwi_dict.items())

for dirnum in set(dwi_dict['bval']):
idxs = { i for k,v in dwi_dict.iteritems() for i in range(0,len(dwi_dict['bval'])) if v[i] == dirnum }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add a comment explaining what this line does?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Line 450: I made a mistake - it's supposed to be
assert all(len(dwi_dict[k]) ==n for k,v in dwi_dict.items())

Line 453: extracts the index values in dwi_dict['bval'] if the value matches "dirnum", which is the number of directions (i.e. 98 or 99). These index values gets used to find the PE directions, dwi files names, etc. in the dictionary.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you could add this comment to the code it would be easier for people in the future to understand this code. Thanks!

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Of course!

PEdirNums = set([dwi_dict['direction'][i] for i in idxs])
for PEdirNum in PEdirNums:
dwis = [ dwi_dict['dwiFile'][i] for i in idxs if dwi_dict['direction'][i] == PEdirNum ]
assert len(dwis) <= 2
dwiname = "Diffusion" + "_dir-" + str(dirnum) + "_" + session + "_corr_" + str(PEdirNum)
if "j" in PEdirNum:
PEdir = 2
elif "i" in PEdirNum:
PEdir = 1
else:
RuntimeError("Phase encoding direction not specified for diffusion data.")
pos = "EMPTY"
neg = "EMPTY"
gdcoeffs = "None"
for dwi in dwis:
if "-" in layout.get_metadata(dwi)["PhaseEncodingDirection"]:
neg = dwi
else:
pos = dwi
try:
echospacing = layout.get_metadata(pos)["EffectiveEchoSpacing"] * 1000
dwi_stage_dict = OrderedDict([("DiffusionPreprocessing", partial(run_diffusion_processsing,
posData=pos,
negData=neg,
path=args.output_dir,
subject="sub-%s" % subject_label,
echospacing=echospacing,
PEdir=PEdir,
gdcoeffs="NONE",
dwiname=dwiname,
n_cpus=args.n_cpus))])
for stage, stage_func in dwi_stage_dict.iteritems():
if stage in args.stages:
stage_func()
except:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please replace this with a more specific catch (narrower in the amount of code it encomapses and types of exceptions it is suppose to catch). Right now it will catch all exceptions making debugging really hard.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Got it!

print("You may have missing diffusion data in the positive phase encoding direction.")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This tra/except statement still encompases a huge chunk of code that could produce "NameError" that does not necessarily mean that there is missing diffusion data in the positive phase encoding direction. This should be made more specific. Furthermore we should exit with a non zero code to indicate that something did not work correctly.