fsl_mrs 15.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.auxiliary 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
    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).')
102
103
    optional.add_argument('--h2o_scale', type=float, default=1.0,
                          help='Additional scaling modifier for external water referencing.')
William Clarke's avatar
William Clarke committed
104
    optional.add_argument('--report', action="store_true",
105
                          help='output html report')
William Clarke's avatar
William Clarke committed
106
    optional.add_argument('--verbose', action="store_true",
107
                          help='spit out verbose info')
William Clarke's avatar
William Clarke committed
108
    optional.add_argument('--phase_correct', action="store_true",
Saad Jbabdi's avatar
Saad Jbabdi committed
109
                          help='do phase correction')
William Clarke's avatar
William Clarke committed
110
    optional.add_argument('--overwrite', action="store_true",
111
                          help='overwrite existing output folder')
William Clarke's avatar
William Clarke committed
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    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
128
129
130

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

132
133
134
135
136
137
    # 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
138
139
140
    import time
    import os
    import shutil
141
    import json
William Clarke's avatar
William Clarke committed
142
    import warnings
William Clarke's avatar
William Clarke committed
143
144
    import matplotlib
    matplotlib.use('agg')
145
146
147
148
149
150
151
    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
        f.write(json.dumps(vars(args)))
175
176
        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 = 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

Saad Jbabdi's avatar
Saad Jbabdi committed
194
    # Instantiate MRS object
William Clarke's avatar
William Clarke committed
195
196
197
198
199
200
201
202
    mrs = FID.mrs(basis=basis,
                  names=names,
                  basis_hdr=basisheader[0],
                  ref_data=H2O)

    if isinstance(mrs, list):
        raise ValueError('fsl_mrs only handles a single FID at a time.'
                         ' Please preprocess data.')
William Clarke's avatar
William Clarke committed
203

204
205
206
207
208
    # 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
209
        conjugated = mrs.check_FID(repair=True)
210
211
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
212
213
                warnings.warn('FID has been checked and conjugated.'
                              ' Please check!', UserWarning)
214
215
216
217
218

    if args.conjbasis is not None:
        if args.conjbasis:
            mrs.conj_Basis()
    else:
William Clarke's avatar
William Clarke committed
219
        conjugated = mrs.check_Basis(repair=True)
220
221
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
222
223
224
                warnings.warn('Basis has been checked and conjugated.'
                              ' Please check!', UserWarning)

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

229
    # Do phase correction
Saad Jbabdi's avatar
Saad Jbabdi committed
230
231
232
    if args.phase_correct:
        if args.verbose:
            print('--->>  Phase correction\n')
William Clarke's avatar
William Clarke committed
233
234
        mrs.FID = misc.phase_correct(mrs, mrs.FID)

235
    # Keep/Ignore metabolites
236
    mrs.keep(args.keep)
William Clarke's avatar
William Clarke committed
237
    mrs.ignore(args.ignore)
238

Saad Jbabdi's avatar
Saad Jbabdi committed
239
240
241
    # Do the fitting here
    if args.verbose:
        print('--->> Start fitting\n\n')
242
        print('    Algorithm = [{}]\n'.format(args.algo))
William Clarke's avatar
William Clarke committed
243
    start = time.time()
Saad Jbabdi's avatar
Saad Jbabdi committed
244

William Clarke's avatar
William Clarke committed
245
    ppmlim = args.ppmlim
246
    if ppmlim is not None:
William Clarke's avatar
William Clarke committed
247
        ppmlim = (ppmlim[0], ppmlim[1])
248
249

    # Parse metabolite groups
William Clarke's avatar
William Clarke committed
250
    metab_groups = misc.parse_metab_groups(mrs, args.metab_groups)
251
252

    # Include Macromolecules? These should have their own metab groups
253
    if args.add_MM:
William Clarke's avatar
William Clarke committed
254
        if args.verbose:
Saad Jbabdi's avatar
Saad Jbabdi committed
255
            print('Adding macromolecules')
William Clarke's avatar
William Clarke committed
256
        nMM = mrs.add_MM_peaks(gamma=10, sigma=20)
William Clarke's avatar
William Clarke committed
257
        G = [i + max(metab_groups) + 1 for i in range(nMM)]
258
259
        metab_groups += G

William Clarke's avatar
William Clarke committed
260
    # Choose fitting lineshape model.
261
    if args.lorentzian:
William Clarke's avatar
William Clarke committed
262
263
264
265
        Fitargs = {'ppmlim': ppmlim,
                   'method': args.algo,
                   'baseline_order': args.baseline_order,
                   'metab_groups': metab_groups,
William Clarke's avatar
William Clarke committed
266
267
                   'model': 'lorentzian',
                   'disable_mh_priors': args.disable_MH_priors}
268
    else:
William Clarke's avatar
William Clarke committed
269
270
271
272
        Fitargs = {'ppmlim': ppmlim,
                   'method': args.algo,
                   'baseline_order': args.baseline_order,
                   'metab_groups': metab_groups,
William Clarke's avatar
William Clarke committed
273
274
                   'model': 'voigt',
                   'disable_mh_priors': args.disable_MH_priors}
275

Saad Jbabdi's avatar
Saad Jbabdi committed
276
277
278
279
280
    if args.verbose:
        print(mrs)
        print('Fitting args:')
        print(Fitargs)

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

283
284
285
    # Quantification
    # Echo time
    if args.TE is not None:
William Clarke's avatar
William Clarke committed
286
        echotime = args.TE * 1E-3
William Clarke's avatar
William Clarke committed
287
288
289
290
291
292
293
    elif 'meta' in basisheader and 'TE' in basisheader['meta']:
        echotime = basisheader['meta']['TE']
        if echotime > 1.0:  # Assume in ms.
            echotime *= 1E-3
        else:
            echotime = None
    elif 'EchoTime' in FID.hdr_ext:
294
        echotime = FID.hdr_ext['EchoTime']
295
296
    else:
        echotime = None
William Clarke's avatar
William Clarke committed
297

298
299
    # Internal and Water quantification if requested
    if (mrs.H2O is None) or (echotime is None):
300
        if echotime is None:
William Clarke's avatar
William Clarke committed
301
302
303
304
            warnings.warn('H2O file provided but could not determine TE:'
                          ' no absolute quantification will be performed.',
                          UserWarning)
        res.calculateConcScaling(mrs, referenceMetab=args.internal_ref)
305
    elif args.tissue_frac is not None:
306
        res.calculateConcScaling(mrs,
William Clarke's avatar
William Clarke committed
307
308
                                 referenceMetab=args.internal_ref,
                                 waterRefFID=mrs.H2O,
309
                                 tissueFractions=args.tissue_frac,
William Clarke's avatar
William Clarke committed
310
                                 TE=echotime,
311
312
                                 verbose=args.verbose,
                                 add_scale=args.h2o_scale)
313
    else:
314
        res.calculateConcScaling(mrs,
William Clarke's avatar
William Clarke committed
315
316
                                 referenceMetab=args.internal_ref,
                                 waterRefFID=mrs.H2O,
317
                                 tissueFractions=None,
William Clarke's avatar
William Clarke committed
318
                                 TE=echotime,
319
320
321
                                 verbose=args.verbose,
                                 add_scale=args.h2o_scale)

322
323
324
    # Combine metabolites.
    if args.combine is not None:
        res.combine(args.combine)
Saad Jbabdi's avatar
Saad Jbabdi committed
325
326
327
328
    stop = time.time()

    # Report on the fitting
    if args.verbose:
William Clarke's avatar
William Clarke committed
329
        duration = stop - start
William Clarke's avatar
William Clarke committed
330
        print(f'    Fitting lasted          : {duration:.3f} secs.\n')
Saad Jbabdi's avatar
Saad Jbabdi committed
331
332
333
    # Save output files
    if args.verbose:
        print('--->> Saving output files to {}\n'.format(args.output))
334

William Clarke's avatar
William Clarke committed
335
336
337
338
339
340
341
342
    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')
343
    if args.algo == 'MH':
William Clarke's avatar
William Clarke committed
344
345
346
347
348
        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
349

Saad Jbabdi's avatar
Saad Jbabdi committed
350
    # Save image of MRS voxel
351
    location_fig = None
352
353
354
355
356
357
    if args.t1 is not None \
            and FID.getXFormCode() > 0:
        fig = plotting.plot_world_orient(args.t1, args.data)
        fig.tight_layout()
        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
358
359

    # Save quick summary figure
William Clarke's avatar
William Clarke committed
360
361
362
363
    report.fitting_summary_fig(mrs, res,
                               filename=os.path.join(args.output,
                                                     'fit_summary.png'))

Saad Jbabdi's avatar
Saad Jbabdi committed
364
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
365
    if args.report:
William Clarke's avatar
William Clarke committed
366
367
368
369
370
371
372
373
374
375
        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
376
377
378
    if args.verbose:
        print('\n\n\nDone.')

William Clarke's avatar
William Clarke committed
379
380
381

def str_or_int_arg(x):
    try:
382
        return int(x)
William Clarke's avatar
William Clarke committed
383
    except ValueError:
384
        return x
385

William Clarke's avatar
William Clarke committed
386

387
def tissue_frac_arg(x):
Saad Jbabdi's avatar
Saad Jbabdi committed
388
    import json
389
390
    try:
        with open(x) as jsonFile:
William Clarke's avatar
William Clarke committed
391
            jsonString = jsonFile.read()
392
        return json.loads(jsonString)
William Clarke's avatar
William Clarke committed
393
    except IOError:
394
395
        return float(x)

William Clarke's avatar
William Clarke committed
396

397
398
399
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
400
401
        if isinstance(values[0], dict):
            setattr(namespace, self.dest, values[0])
402
        else:
William Clarke's avatar
William Clarke committed
403
404
405
406
            setattr(namespace, self.dest,
                    {'WM': values[0], 'GM': values[1], 'CSF': values[2]})


Saad Jbabdi's avatar
Saad Jbabdi committed
407
408
if __name__ == '__main__':
    main()