fsl_mrs 16.5 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)')
William Clarke's avatar
William Clarke committed
80
81
82
83
84
85
    fitting_args.add_argument('--ind_scale', default=None, type=str,
                              nargs='+',
                              help='List of basis spectra to scale'
                                   ' independently of other basis spectra.')
    fitting_args.add_argument('--disable_MH_priors', action="store_true",
                              help="Disable MH priors.")
86

William Clarke's avatar
William Clarke committed
87
88
    # ADDITIONAL OPTIONAL ARGUMENTS
    optional.add_argument('--t1', type=str, default=None, metavar='IMAGE',
Saad Jbabdi's avatar
Saad Jbabdi committed
89
                          help='structural image (for report)')
William Clarke's avatar
William Clarke committed
90
    optional.add_argument('--TE', type=float, default=None, metavar='TE',
91
                          help='Echo time for relaxation correction (ms)')
William Clarke's avatar
William Clarke committed
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
    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",
107
                          help='output html report')
William Clarke's avatar
William Clarke committed
108
    optional.add_argument('--verbose', action="store_true",
109
                          help='spit out verbose info')
William Clarke's avatar
William Clarke committed
110
    optional.add_argument('--phase_correct', action="store_true",
Saad Jbabdi's avatar
Saad Jbabdi committed
111
                          help='do phase correction')
William Clarke's avatar
William Clarke committed
112
    optional.add_argument('--overwrite', action="store_true",
113
                          help='overwrite existing output folder')
William Clarke's avatar
William Clarke committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    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
130
131
132

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

134
135
136
137
138
139
    # 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
140
141
142
143
    import time
    import os
    import shutil
    import warnings
144
145
146
147
148
149
150
151
    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
152
153
    if not args.verbose:
        warnings.filterwarnings("ignore")
154

Saad Jbabdi's avatar
Saad Jbabdi committed
155
156
157
158
    # Check if output folder exists
    overwrite = args.overwrite
    if os.path.exists(args.output):
        if not overwrite:
William Clarke's avatar
William Clarke committed
159
160
            print(f"Folder '{args.output}' exists."
                  " Are you sure you want to delete it? [Y,N]")
Saad Jbabdi's avatar
Saad Jbabdi committed
161
162
163
164
165
166
167
            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
168
            os.makedirs(args.output, exist_ok=True)
Saad Jbabdi's avatar
Saad Jbabdi committed
169
    else:
William Clarke's avatar
William Clarke committed
170
        os.makedirs(args.output, exist_ok=True)
171
172

    # Save chosen arguments
William Clarke's avatar
William Clarke committed
173
    with open(os.path.join(args.output, "options.txt"), "w") as f:
174
175
176
        f.write(str(args))
        f.write("\n--------\n")
        f.write(p.format_values())
William Clarke's avatar
William Clarke committed
177
178

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

Saad Jbabdi's avatar
Saad Jbabdi committed
180
    # Read data/h2o/basis
William Clarke's avatar
William Clarke committed
181
    if args.verbose:
Saad Jbabdi's avatar
Saad Jbabdi committed
182
        print('--->> Read input data and basis\n')
William Clarke's avatar
William Clarke committed
183
184
        print(f'  {args.data}')
        print(f'  {args.basis}\n')
Saad Jbabdi's avatar
Saad Jbabdi committed
185

William Clarke's avatar
William Clarke committed
186
    FID, dataheader = mrs_io.read_FID(args.data)
187
    basis, names, basisheader = mrs_io.read_basis(args.basis)
William Clarke's avatar
William Clarke committed
188

189
    if args.h2o is not None:
William Clarke's avatar
William Clarke committed
190
        H2O, _ = mrs_io.read_FID(args.h2o)
191
192
    else:
        H2O = None
William Clarke's avatar
William Clarke committed
193

194
195
196
197
198
    # 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
199
        if args.verbose:
William Clarke's avatar
William Clarke committed
200
201
            print(f'     Detected central frequency'
                  f' in header info cf = {cf:0.6f} MHz')
202
    else:
William Clarke's avatar
William Clarke committed
203
204
        raise(Exception('Cannot determine central frequency.'
                        'Please either set it or include it in data header'))
205
206
207
208
209

    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
210
        if args.verbose:
211
            print(f'      Detected bandwidth in header info bw = {bw:0.1f} Hz')
212
    else:
William Clarke's avatar
William Clarke committed
213
214
        raise(Exception('Cannot determine bandwidth.'
                        'Please either set it or include it in data header'))
Saad Jbabdi's avatar
Saad Jbabdi committed
215
216
217

    # Fix case where basis file contains no header info (e.g. .RAW)
    if basisheader is None:
William Clarke's avatar
William Clarke committed
218
219
220
        basisheader = {'bandwidth': bw,
                       'dwelltime': 1/bw,
                       'centralFrequency': cf}
Saad Jbabdi's avatar
Saad Jbabdi committed
221
222
    else:
        basisheader = basisheader[0]
William Clarke's avatar
William Clarke committed
223

Saad Jbabdi's avatar
Saad Jbabdi committed
224
    # Instantiate MRS object
William Clarke's avatar
William Clarke committed
225
226
227
228
229
    MRSargs = {'FID': FID, 'basis': basis,
               'basis_hdr': basisheader, 'names': names,
               'H2O': H2O, 'cf': cf, 'bw': bw}
    mrs = MRS(**MRSargs)

230
231
232
233
234
    # 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
235
        conjugated = mrs.check_FID(repair=True)
236
237
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
238
239
                warnings.warn('FID has been checked and conjugated.'
                              ' Please check!', UserWarning)
240
241
242
243
244

    if args.conjbasis is not None:
        if args.conjbasis:
            mrs.conj_Basis()
    else:
William Clarke's avatar
William Clarke committed
245
        conjugated = mrs.check_Basis(repair=True)
246
247
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
248
249
250
                warnings.warn('Basis has been checked and conjugated.'
                              ' Please check!', UserWarning)

251
252
    # Rescale FID, H2O and basis to have nice range
    if not args.no_rescale:
William Clarke's avatar
William Clarke committed
253
        mrs.rescaleForFitting(ind_scaling=args.ind_scale)
Saad Jbabdi's avatar
Saad Jbabdi committed
254

255
    # Do phase correction
Saad Jbabdi's avatar
Saad Jbabdi committed
256
257
258
    if args.phase_correct:
        if args.verbose:
            print('--->>  Phase correction\n')
William Clarke's avatar
William Clarke committed
259
260
        mrs.FID = misc.phase_correct(mrs, mrs.FID)

261
    # Keep/Ignore metabolites
262
    mrs.keep(args.keep)
William Clarke's avatar
William Clarke committed
263
    mrs.ignore(args.ignore)
264

Saad Jbabdi's avatar
Saad Jbabdi committed
265
266
267
    # Do the fitting here
    if args.verbose:
        print('--->> Start fitting\n\n')
268
        print('    Algorithm = [{}]\n'.format(args.algo))
Saad Jbabdi's avatar
Saad Jbabdi committed
269
270
        start = time.time()

William Clarke's avatar
William Clarke committed
271
    ppmlim = args.ppmlim
272
    if ppmlim is not None:
William Clarke's avatar
William Clarke committed
273
        ppmlim = (ppmlim[0], ppmlim[1])
274
275

    # Parse metabolite groups
William Clarke's avatar
William Clarke committed
276
    metab_groups = misc.parse_metab_groups(mrs, args.metab_groups)
277
278

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

William Clarke's avatar
William Clarke committed
286
    # Choose fitting lineshape model.
287
    if args.lorentzian:
William Clarke's avatar
William Clarke committed
288
289
290
291
        Fitargs = {'ppmlim': ppmlim,
                   'method': args.algo,
                   'baseline_order': args.baseline_order,
                   'metab_groups': metab_groups,
William Clarke's avatar
William Clarke committed
292
293
                   'model': 'lorentzian',
                   'disable_mh_priors': args.disable_MH_priors}
294
    else:
William Clarke's avatar
William Clarke committed
295
296
297
298
        Fitargs = {'ppmlim': ppmlim,
                   'method': args.algo,
                   'baseline_order': args.baseline_order,
                   'metab_groups': metab_groups,
William Clarke's avatar
William Clarke committed
299
300
                   'model': 'voigt',
                   'disable_mh_priors': args.disable_MH_priors}
301

Saad Jbabdi's avatar
Saad Jbabdi committed
302
303
304
305
306
    if args.verbose:
        print(mrs)
        print('Fitting args:')
        print(Fitargs)

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

309
310
311
312
    # Quantification
    # Echo time
    if args.TE is not None:
        echotime = args.TE*1E-3
313
314
315
    elif 'meta' in basisheader:
        if 'TE' in basisheader['meta']:
            echotime = basisheader['meta']['TE']
William Clarke's avatar
William Clarke committed
316
            if echotime > 1.0:  # Assume in ms.
317
318
319
                echotime *= 1E-3
            else:
                echotime = None
320
321
    elif 'TE' in dataheader:
        echotime = dataheader['TE']
322
323
    else:
        echotime = None
William Clarke's avatar
William Clarke committed
324

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

    # Report on the fitting
    if args.verbose:
William Clarke's avatar
William Clarke committed
353
354
        duration = stop-start
        print(f'    Fitting lasted          : {duration:.3f} secs.\n')
Saad Jbabdi's avatar
Saad Jbabdi committed
355
356
357
    # Save output files
    if args.verbose:
        print('--->> Saving output files to {}\n'.format(args.output))
358

William Clarke's avatar
William Clarke committed
359
360
361
362
363
364
365
366
    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')
367
    if args.algo == 'MH':
William Clarke's avatar
William Clarke committed
368
369
370
371
372
        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
373

Saad Jbabdi's avatar
Saad Jbabdi committed
374
    # Save image of MRS voxel
375
    location_fig = None
Saad Jbabdi's avatar
Saad Jbabdi committed
376
377
    if args.t1 is not None:
        datatype = mrs_io.check_datatype(args.data)
William Clarke's avatar
William Clarke committed
378
379
        if datatype == 'NIFTI':
            fig = plotting.plot_world_orient(args.t1, args.data)
380
            fig.tight_layout()
William Clarke's avatar
William Clarke committed
381
382
            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
383
384

    # Save quick summary figure
William Clarke's avatar
William Clarke committed
385
386
387
388
    report.fitting_summary_fig(mrs, res,
                               filename=os.path.join(args.output,
                                                     'fit_summary.png'))

Saad Jbabdi's avatar
Saad Jbabdi committed
389
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
390
    if args.report:
William Clarke's avatar
William Clarke committed
391
392
393
394
395
396
397
398
399
400
        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
401
402
403
    if args.verbose:
        print('\n\n\nDone.')

William Clarke's avatar
William Clarke committed
404
405
406

def str_or_int_arg(x):
    try:
407
        return int(x)
William Clarke's avatar
William Clarke committed
408
    except ValueError:
409
        return x
410

William Clarke's avatar
William Clarke committed
411

412
def tissue_frac_arg(x):
Saad Jbabdi's avatar
Saad Jbabdi committed
413
    import json
414
415
    try:
        with open(x) as jsonFile:
William Clarke's avatar
William Clarke committed
416
            jsonString = jsonFile.read()
417
        return json.loads(jsonString)
William Clarke's avatar
William Clarke committed
418
    except IOError:
419
420
        return float(x)

William Clarke's avatar
William Clarke committed
421

422
423
424
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
425
426
        if isinstance(values[0], dict):
            setattr(namespace, self.dest, values[0])
427
        else:
William Clarke's avatar
William Clarke committed
428
429
430
431
            setattr(namespace, self.dest,
                    {'WM': values[0], 'GM': values[1], 'CSF': values[2]})


Saad Jbabdi's avatar
Saad Jbabdi committed
432
433
if __name__ == '__main__':
    main()