fsl_mrs 16 KB
Newer Older
Saad Jbabdi's avatar
Saad Jbabdi committed
1
2
3
4
5
#!/usr/bin/env python

# fsl_mrs - wrapper script for MRS fitting
#
# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
6
#         William Carke <william.clarke@ndcn.ox.ac.uk>
Saad Jbabdi's avatar
Saad Jbabdi committed
7
#
William Clarke's avatar
William Clarke committed
8
# Copyright (C) 2019 University of Oxford
Saad Jbabdi's avatar
Saad Jbabdi committed
9
10
# SHBASECOPYRIGHT

11
# Quick imports
12
from fsl_mrs.aux import configargparse
Saad Jbabdi's avatar
Saad Jbabdi committed
13

14
15
from fsl_mrs import __version__
from fsl_mrs.utils.splash import splash
Saad Jbabdi's avatar
Saad Jbabdi committed
16

Saad Jbabdi's avatar
Saad Jbabdi committed
17

18
# NOTE!!!! THERE ARE MORE IMPORTS IN THE CODE BELOW (AFTER ARGPARSING)
Saad Jbabdi's avatar
Saad Jbabdi committed
19
20
21

def main():
    # Parse command-line arguments
William Clarke's avatar
William Clarke committed
22
23
24
    p = configargparse.ArgParser(
        add_config_file_help=False,
        description="FSL Magnetic Resonance Spectroscopy Wrapper Script")
25
26
27
28

    # utility for hiding certain arguments
    def hide_args(arglist):
        for action in arglist:
William Clarke's avatar
William Clarke committed
29
30
31
            action.help = p.SUPPRESS

    p.add_argument('-v', '--version', action='version', version=__version__)
32

William Clarke's avatar
William Clarke committed
33
    required = p.add_argument_group('required arguments')
Saad Jbabdi's avatar
Saad Jbabdi committed
34
    fitting_args = p.add_argument_group('fitting options')
William Clarke's avatar
William Clarke committed
35
    optional = p.add_argument_group('additional options')
Saad Jbabdi's avatar
Saad Jbabdi committed
36
37

    # REQUIRED ARGUMENTS
38
    required.add_argument('--data',
William Clarke's avatar
William Clarke committed
39
                          required=True, type=str, metavar='<str>',
Saad Jbabdi's avatar
Saad Jbabdi committed
40
                          help='input FID file')
41
    required.add_argument('--basis',
William Clarke's avatar
William Clarke committed
42
43
44
                          required=True, type=str, metavar='<str>',
                          help='.BASIS file or folder containing basis spectra'
                               '(will read all files within)')
45
    required.add_argument('--output',
William Clarke's avatar
William Clarke committed
46
                          required=True, type=str, metavar='<str>',
Saad Jbabdi's avatar
Saad Jbabdi committed
47
                          help='output folder')
William Clarke's avatar
William Clarke committed
48

Saad Jbabdi's avatar
Saad Jbabdi committed
49
    # FITTING ARGUMENTS
William Clarke's avatar
William Clarke committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    fitting_args.add_argument('--algo', default='Newton', type=str,
                              help='algorithm [Newton (fast, default)'
                                   ' or MH (slow)]')
    fitting_args.add_argument('--ignore', type=str, nargs='+',
                              metavar='METAB',
                              help='ignore certain metabolites [repeatable]')
    fitting_args.add_argument('--keep', type=str, nargs='+', metavar='METAB',
                              help='only keep these metabolites')
    fitting_args.add_argument('--combine', type=str, nargs='+',
                              action='append', metavar='METAB',
                              help='combine certain metabolites [repeatable]')
    fitting_args.add_argument('--ppmlim', default=(.2, 4.2), type=float,
                              nargs=2, metavar=('LOW', 'HIGH'),
                              help='limit the fit to a freq range'
                                   ' (default=(.2,4.2))')
    fitting_args.add_argument('--h2o', default=None, type=str, metavar='H2O',
                              help='input .H2O file for quantification')
    fitting_args.add_argument('--baseline_order', default=2, type=int,
                              metavar=('ORDER'),
                              help='order of baseline polynomial'
                                   ' (default=2,-1 disables)')
    fitting_args.add_argument('--metab_groups', default=0, nargs='+',
                              type=str_or_int_arg,
                              help='metabolite groups: list of groups'
                                   ' or list of names for indept groups.')
    fitting_args.add_argument('--add_MM', action="store_true",
76
                              help="include default macromolecule peaks")
William Clarke's avatar
William Clarke committed
77
78
79
    fitting_args.add_argument('--lorentzian', action="store_true",
                              help='Enable purely lorentzian broadening'
                                   ' (default is Voigt)')
80

William Clarke's avatar
William Clarke committed
81
82
    # ADDITIONAL OPTIONAL ARGUMENTS
    optional.add_argument('--t1', type=str, default=None, metavar='IMAGE',
Saad Jbabdi's avatar
Saad Jbabdi committed
83
                          help='structural image (for report)')
William Clarke's avatar
William Clarke committed
84
    optional.add_argument('--TE', type=float, default=None, metavar='TE',
85
                          help='Echo time for relaxation correction (ms)')
William Clarke's avatar
William Clarke committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    optional.add_argument('--tissue_frac', type=tissue_frac_arg,
                          action=TissueFracAction, nargs='+',
                          default=None, metavar='WM GM CSF OR json',
                          help='Fractional tissue volumes for WM, GM, CSF'
                               ' or json segmentation file.'
                               ' Defaults to pure water scaling.')
    optional.add_argument('--internal_ref', type=str, default=['Cr', 'PCr'],
                          nargs='+',
                          help='Metabolite(s) used as an internal reference.'
                               ' Defaults to tCr (Cr+PCr).')
    optional.add_argument('--central_frequency', default=None, type=float,
                          help='central frequency in Hz')
    optional.add_argument('--dwell_time', default=None, type=float,
                          help='dwell time in seconds')
    optional.add_argument('--report', action="store_true",
101
                          help='output html report')
William Clarke's avatar
William Clarke committed
102
    optional.add_argument('--verbose', action="store_true",
103
                          help='spit out verbose info')
William Clarke's avatar
William Clarke committed
104
    optional.add_argument('--phase_correct', action="store_true",
Saad Jbabdi's avatar
Saad Jbabdi committed
105
                          help='do phase correction')
William Clarke's avatar
William Clarke committed
106
    optional.add_argument('--overwrite', action="store_true",
107
                          help='overwrite existing output folder')
William Clarke's avatar
William Clarke committed
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    optional.add_argument('--conj_fid', dest='conjfid', action="store_true",
                          help='Force conjugation of FID')
    optional.add_argument('--no_conj_fid', dest='conjfid',
                          action="store_false",
                          help='Forbid automatic conjugation of FID')
    optional.add_argument('--conj_basis', dest='conjbasis',
                          action="store_true",
                          help='Force conjugation of basis')
    optional.add_argument('--no_conj_basis', dest='conjbasis',
                          action="store_false",
                          help='Forbid automatic conjugation of basis')
    optional.set_defaults(conjfid=None, conjbasis=None)
    optional.add_argument('--no_rescale', action="store_true",
                          help='Forbid rescaling of FID/basis/H2O.')
    optional.add('--config', required=False, is_config_file=True,
                 help='configuration file')
Saad Jbabdi's avatar
Saad Jbabdi committed
124
125
126

    # Parse command-line arguments
    args = p.parse_args()
William Clarke's avatar
William Clarke committed
127

128
129
130
131
132
133
    # Output kickass splash screen
    if args.verbose:
        splash(logo='mrs')

    # ######################################################
    # DO THE IMPORTS AFTER PARSING TO SPEED UP HELP DISPLAY
William Clarke's avatar
William Clarke committed
134
135
136
137
    import time
    import os
    import shutil
    import warnings
138
139
140
141
142
143
144
145
    from fsl_mrs.core import MRS
    from fsl_mrs.utils import mrs_io
    from fsl_mrs.utils import report
    from fsl_mrs.utils import fitting
    from fsl_mrs.utils import plotting
    from fsl_mrs.utils import misc
    import datetime
    # ######################################################
Saad Jbabdi's avatar
Saad Jbabdi committed
146
147
    if not args.verbose:
        warnings.filterwarnings("ignore")
148

Saad Jbabdi's avatar
Saad Jbabdi committed
149
150
151
152
    # Check if output folder exists
    overwrite = args.overwrite
    if os.path.exists(args.output):
        if not overwrite:
William Clarke's avatar
William Clarke committed
153
154
            print(f"Folder '{args.output}' exists."
                  " Are you sure you want to delete it? [Y,N]")
Saad Jbabdi's avatar
Saad Jbabdi committed
155
156
157
158
159
160
161
            response = input()
            overwrite = response.upper() == "Y"
        if not overwrite:
            print('Early stopping...')
            exit()
        else:
            shutil.rmtree(args.output)
William Clarke's avatar
William Clarke committed
162
            os.makedirs(args.output, exist_ok=True)
Saad Jbabdi's avatar
Saad Jbabdi committed
163
    else:
William Clarke's avatar
William Clarke committed
164
        os.makedirs(args.output, exist_ok=True)
165
166

    # Save chosen arguments
William Clarke's avatar
William Clarke committed
167
    with open(os.path.join(args.output, "options.txt"), "w") as f:
168
169
170
        f.write(str(args))
        f.write("\n--------\n")
        f.write(p.format_values())
William Clarke's avatar
William Clarke committed
171
172

    # Do the work
Saad Jbabdi's avatar
Saad Jbabdi committed
173

Saad Jbabdi's avatar
Saad Jbabdi committed
174
    # Read data/h2o/basis
William Clarke's avatar
William Clarke committed
175
    if args.verbose:
Saad Jbabdi's avatar
Saad Jbabdi committed
176
        print('--->> Read input data and basis\n')
William Clarke's avatar
William Clarke committed
177
178
        print(f'  {args.data}')
        print(f'  {args.basis}\n')
Saad Jbabdi's avatar
Saad Jbabdi committed
179

William Clarke's avatar
William Clarke committed
180
    FID, dataheader = mrs_io.read_FID(args.data)
181
    basis, names, basisheader = mrs_io.read_basis(args.basis)
William Clarke's avatar
William Clarke committed
182

183
    if args.h2o is not None:
William Clarke's avatar
William Clarke committed
184
        H2O, _ = mrs_io.read_FID(args.h2o)
185
186
    else:
        H2O = None
William Clarke's avatar
William Clarke committed
187

188
189
190
191
192
    # Collect useful info
    if args.central_frequency is not None:
        cf = args.central_frequency
    elif dataheader['centralFrequency'] is not None:
        cf = dataheader['centralFrequency']
Saad Jbabdi's avatar
Saad Jbabdi committed
193
        if args.verbose:
William Clarke's avatar
William Clarke committed
194
195
            print(f'     Detected central frequency'
                  f' in header info cf = {cf:0.6f} MHz')
196
    else:
William Clarke's avatar
William Clarke committed
197
198
        raise(Exception('Cannot determine central frequency.'
                        'Please either set it or include it in data header'))
199
200
201
202
203

    if args.dwell_time is not None:
        bw = 1/args.dwell_time
    elif dataheader['bandwidth'] is not None:
        bw = dataheader['bandwidth']
Saad Jbabdi's avatar
Saad Jbabdi committed
204
        if args.verbose:
205
            print(f'      Detected bandwidth in header info bw = {bw:0.1f} Hz')
206
    else:
William Clarke's avatar
William Clarke committed
207
208
        raise(Exception('Cannot determine bandwidth.'
                        'Please either set it or include it in data header'))
Saad Jbabdi's avatar
Saad Jbabdi committed
209
210
211

    # Fix case where basis file contains no header info (e.g. .RAW)
    if basisheader is None:
William Clarke's avatar
William Clarke committed
212
213
214
        basisheader = {'bandwidth': bw,
                       'dwelltime': 1/bw,
                       'centralFrequency': cf}
Saad Jbabdi's avatar
Saad Jbabdi committed
215
216
    else:
        basisheader = basisheader[0]
William Clarke's avatar
William Clarke committed
217

Saad Jbabdi's avatar
Saad Jbabdi committed
218
    # Instantiate MRS object
William Clarke's avatar
William Clarke committed
219
220
221
222
223
    MRSargs = {'FID': FID, 'basis': basis,
               'basis_hdr': basisheader, 'names': names,
               'H2O': H2O, 'cf': cf, 'bw': bw}
    mrs = MRS(**MRSargs)

224
225
226
227
228
    # Check the FID and basis / conjugate
    if args.conjfid is not None:
        if args.conjfid:
            mrs.conj_FID()
    else:
William Clarke's avatar
William Clarke committed
229
        conjugated = mrs.check_FID(repair=True)
230
231
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
232
233
                warnings.warn('FID has been checked and conjugated.'
                              ' Please check!', UserWarning)
234
235
236
237
238

    if args.conjbasis is not None:
        if args.conjbasis:
            mrs.conj_Basis()
    else:
William Clarke's avatar
William Clarke committed
239
        conjugated = mrs.check_Basis(repair=True)
240
241
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
242
243
244
                warnings.warn('Basis has been checked and conjugated.'
                              ' Please check!', UserWarning)

245
246
247
    # Rescale FID, H2O and basis to have nice range
    if not args.no_rescale:
        mrs.rescaleForFitting()
Saad Jbabdi's avatar
Saad Jbabdi committed
248

249
    # Do phase correction
Saad Jbabdi's avatar
Saad Jbabdi committed
250
251
252
    if args.phase_correct:
        if args.verbose:
            print('--->>  Phase correction\n')
William Clarke's avatar
William Clarke committed
253
        mrs.FID = misc.phase_correct(mrs, mrs.FID)
254
        mrs.Spec = misc.FIDToSpec(mrs.FID)
William Clarke's avatar
William Clarke committed
255

256
    # Keep/Ignore metabolites
257
    mrs.keep(args.keep)
William Clarke's avatar
William Clarke committed
258
    mrs.ignore(args.ignore)
259

Saad Jbabdi's avatar
Saad Jbabdi committed
260
261
262
    # Do the fitting here
    if args.verbose:
        print('--->> Start fitting\n\n')
263
        print('    Algorithm = [{}]\n'.format(args.algo))
Saad Jbabdi's avatar
Saad Jbabdi committed
264
265
        start = time.time()

William Clarke's avatar
William Clarke committed
266
    ppmlim = args.ppmlim
267
    if ppmlim is not None:
William Clarke's avatar
William Clarke committed
268
        ppmlim = (ppmlim[0], ppmlim[1])
269
270

    # Parse metabolite groups
William Clarke's avatar
William Clarke committed
271
    metab_groups = misc.parse_metab_groups(mrs, args.metab_groups)
272
273

    # Include Macromolecules? These should have their own metab groups
274
    if args.add_MM:
Saad Jbabdi's avatar
Saad Jbabdi committed
275
276
        if not args.verbose:
            print('Adding macromolecules')
William Clarke's avatar
William Clarke committed
277
278
        nMM = mrs.add_MM_peaks(gamma=10, sigma=20)
        G = [i+max(metab_groups)+1 for i in range(nMM)]
279
280
        metab_groups += G

281
    if args.lorentzian:
William Clarke's avatar
William Clarke committed
282
283
284
285
286
        Fitargs = {'ppmlim': ppmlim,
                   'method': args.algo,
                   'baseline_order': args.baseline_order,
                   'metab_groups': metab_groups,
                   'model': 'lorentzian'}
287
    else:
William Clarke's avatar
William Clarke committed
288
289
290
291
292
        Fitargs = {'ppmlim': ppmlim,
                   'method': args.algo,
                   'baseline_order': args.baseline_order,
                   'metab_groups': metab_groups,
                   'model': 'voigt'}
293

Saad Jbabdi's avatar
Saad Jbabdi committed
294
295
296
297
298
    if args.verbose:
        print(mrs)
        print('Fitting args:')
        print(Fitargs)

William Clarke's avatar
William Clarke committed
299
    res = fitting.fit_FSLModel(mrs, **Fitargs)
300

301
302
303
304
    # Quantification
    # Echo time
    if args.TE is not None:
        echotime = args.TE*1E-3
305
306
307
    elif 'meta' in basisheader:
        if 'TE' in basisheader['meta']:
            echotime = basisheader['meta']['TE']
William Clarke's avatar
William Clarke committed
308
            if echotime > 1.0:  # Assume in ms.
309
310
311
                echotime *= 1E-3
            else:
                echotime = None
312
313
    elif 'TE' in dataheader:
        echotime = dataheader['TE']
314
315
    else:
        echotime = None
William Clarke's avatar
William Clarke committed
316

317
318
    # Internal and Water quantification if requested
    if (mrs.H2O is None) or (echotime is None):
319
        if echotime is None:
William Clarke's avatar
William Clarke committed
320
321
322
323
            warnings.warn('H2O file provided but could not determine TE:'
                          ' no absolute quantification will be performed.',
                          UserWarning)
        res.calculateConcScaling(mrs, referenceMetab=args.internal_ref)
324
    elif args.tissue_frac is None:
325
        res.calculateConcScaling(mrs,
William Clarke's avatar
William Clarke committed
326
327
328
329
330
                                 referenceMetab=args.internal_ref,
                                 waterRefFID=mrs.H2O,
                                 tissueFractions=None,
                                 TE=echotime,
                                 verbose=args.verbose)
331
    else:
332
        res.calculateConcScaling(mrs,
William Clarke's avatar
William Clarke committed
333
334
335
336
337
                                 referenceMetab=args.internal_ref,
                                 waterRefFID=mrs.H2O,
                                 tissueFractions=args.tissue_frac,
                                 TE=echotime,
                                 verbose=args.verbose)
338
339
340
    # Combine metabolites.
    if args.combine is not None:
        res.combine(args.combine)
Saad Jbabdi's avatar
Saad Jbabdi committed
341
342
343
344
    stop = time.time()

    # Report on the fitting
    if args.verbose:
William Clarke's avatar
William Clarke committed
345
346
        duration = stop-start
        print(f'    Fitting lasted          : {duration:.3f} secs.\n')
Saad Jbabdi's avatar
Saad Jbabdi committed
347
348
349
    # Save output files
    if args.verbose:
        print('--->> Saving output files to {}\n'.format(args.output))
350

William Clarke's avatar
William Clarke committed
351
352
353
354
355
356
357
358
    res.to_file(filename=os.path.join(args.output, 'summary.csv'),
                what='summary')
    res.to_file(filename=os.path.join(args.output, 'concentrations.csv'),
                what='concentrations')
    res.to_file(filename=os.path.join(args.output, 'qc.csv'),
                what='qc')
    res.to_file(filename=os.path.join(args.output, 'all_parameters.csv'),
                what='parameters')
359
    if args.algo == 'MH':
William Clarke's avatar
William Clarke committed
360
361
362
363
364
        res.to_file(filename=os.path.join(
                    args.output, 'concentration_samples.csv'),
                    what='concentrations-mh')
        res.to_file(filename=os.path.join(args.output, 'all_samples.csv'),
                    what='parameters-mh')
Saad Jbabdi's avatar
Saad Jbabdi committed
365

Saad Jbabdi's avatar
Saad Jbabdi committed
366
    # Save image of MRS voxel
367
    location_fig = None
Saad Jbabdi's avatar
Saad Jbabdi committed
368
369
    if args.t1 is not None:
        datatype = mrs_io.check_datatype(args.data)
William Clarke's avatar
William Clarke committed
370
371
        if datatype == 'NIFTI':
            fig = plotting.plot_world_orient(args.t1, args.data)
372
            fig.tight_layout()
William Clarke's avatar
William Clarke committed
373
374
            location_fig = os.path.join(args.output, 'voxel_location.png')
            fig.savefig(location_fig, bbox_inches='tight', facecolor='k')
Saad Jbabdi's avatar
Saad Jbabdi committed
375
376

    # Save quick summary figure
William Clarke's avatar
William Clarke committed
377
378
379
380
    report.fitting_summary_fig(mrs, res,
                               filename=os.path.join(args.output,
                                                     'fit_summary.png'))

Saad Jbabdi's avatar
Saad Jbabdi committed
381
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
382
    if args.report:
William Clarke's avatar
William Clarke committed
383
384
385
386
387
388
389
390
391
392
        report.create_report(
            mrs,
            res,
            filename=os.path.join(args.output, 'report.html'),
            fidfile=args.data,
            basisfile=args.basis,
            h2ofile=args.h2o,
            date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"),
            location_fig=location_fig)

Saad Jbabdi's avatar
Saad Jbabdi committed
393
394
395
    if args.verbose:
        print('\n\n\nDone.')

William Clarke's avatar
William Clarke committed
396
397
398

def str_or_int_arg(x):
    try:
399
        return int(x)
William Clarke's avatar
William Clarke committed
400
    except TypeError:
401
        return x
402

William Clarke's avatar
William Clarke committed
403

404
def tissue_frac_arg(x):
Saad Jbabdi's avatar
Saad Jbabdi committed
405
    import json
406
407
    try:
        with open(x) as jsonFile:
William Clarke's avatar
William Clarke committed
408
            jsonString = jsonFile.read()
409
        return json.loads(jsonString)
William Clarke's avatar
William Clarke committed
410
    except IOError:
411
412
        return float(x)

William Clarke's avatar
William Clarke committed
413

414
415
416
class TissueFracAction(configargparse.Action):
    """Sort out tissue fraction types. Should return dict"""
    def __call__(self, parser, namespace, values, option_string=None):
William Clarke's avatar
William Clarke committed
417
418
        if isinstance(values[0], dict):
            setattr(namespace, self.dest, values[0])
419
        else:
William Clarke's avatar
William Clarke committed
420
421
422
423
            setattr(namespace, self.dest,
                    {'WM': values[0], 'GM': values[1], 'CSF': values[2]})


Saad Jbabdi's avatar
Saad Jbabdi committed
424
425
if __name__ == '__main__':
    main()