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.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
    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'
70
                                   ' (default=2, -1 disables)')
William Clarke's avatar
William Clarke committed
71
72
73
74
75
    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)')
92
93
    optional.add_argument('--TR', type=float, default=None, metavar='TR',
                          help='Repetition time for relaxation correction (s)')
William Clarke's avatar
William Clarke committed
94
95
96
97
98
99
100
101
102
103
    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).')
104
105
106
107
108
109
    optional.add_argument('--wref_metabolite', type=str, default=None,
                          help='Metabolite(s) used as an the reference for water scaling.'
                               ' Uses internal defaults otherwise.')
    optional.add_argument('--ref_protons', type=int, default=None,
                          help='Number of protons that reference metabolite is equivalent to.'
                               ' No effect without setting --wref_metabolite.')
110
    optional.add_argument('--ref_int_limits', type=float, default=None, nargs=2,
111
112
                          help='Reference spectrum integration limits (low, high).'
                               ' No effect without setting --wref_metabolite.')
113
114
    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
115
    optional.add_argument('--report', action="store_true",
116
                          help='output html report')
William Clarke's avatar
William Clarke committed
117
    optional.add_argument('--verbose', action="store_true",
118
                          help='spit out verbose info')
William Clarke's avatar
William Clarke committed
119
    optional.add_argument('--overwrite', action="store_true",
120
                          help='overwrite existing output folder')
William Clarke's avatar
William Clarke committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    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
137
138
139

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

141
142
143
144
145
146
    # 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
147
148
149
    import time
    import os
    import shutil
150
    import json
William Clarke's avatar
William Clarke committed
151
    import warnings
William Clarke's avatar
William Clarke committed
152
153
    import matplotlib
    matplotlib.use('agg')
154
155
156
157
158
    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
159
    from fsl_mrs.utils import quantify
160
161
    import datetime
    # ######################################################
Saad Jbabdi's avatar
Saad Jbabdi committed
162
163
    if not args.verbose:
        warnings.filterwarnings("ignore")
164

Saad Jbabdi's avatar
Saad Jbabdi committed
165
166
167
168
    # Check if output folder exists
    overwrite = args.overwrite
    if os.path.exists(args.output):
        if not overwrite:
William Clarke's avatar
William Clarke committed
169
170
            print(f"Folder '{args.output}' exists."
                  " Are you sure you want to delete it? [Y,N]")
Saad Jbabdi's avatar
Saad Jbabdi committed
171
172
173
174
175
176
177
            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
178
            os.makedirs(args.output, exist_ok=True)
Saad Jbabdi's avatar
Saad Jbabdi committed
179
    else:
William Clarke's avatar
William Clarke committed
180
        os.makedirs(args.output, exist_ok=True)
181
182

    # Save chosen arguments
William Clarke's avatar
William Clarke committed
183
    with open(os.path.join(args.output, "options.txt"), "w") as f:
184
        f.write(json.dumps(vars(args)))
185
186
        f.write("\n--------\n")
        f.write(p.format_values())
William Clarke's avatar
William Clarke committed
187
188

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

Saad Jbabdi's avatar
Saad Jbabdi committed
190
    # Read data/h2o/basis
William Clarke's avatar
William Clarke committed
191
    if args.verbose:
Saad Jbabdi's avatar
Saad Jbabdi committed
192
        print('--->> Read input data and basis\n')
William Clarke's avatar
William Clarke committed
193
194
        print(f'  {args.data}')
        print(f'  {args.basis}\n')
Saad Jbabdi's avatar
Saad Jbabdi committed
195

William Clarke's avatar
William Clarke committed
196
    FID = mrs_io.read_FID(args.data)
197
    basis, names, basisheader = mrs_io.read_basis(args.basis)
William Clarke's avatar
William Clarke committed
198

199
    if args.h2o is not None:
William Clarke's avatar
William Clarke committed
200
        H2O = mrs_io.read_FID(args.h2o)
201
202
    else:
        H2O = None
William Clarke's avatar
William Clarke committed
203

Saad Jbabdi's avatar
Saad Jbabdi committed
204
    # Instantiate MRS object
William Clarke's avatar
William Clarke committed
205
206
207
208
209
210
211
212
    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
213

214
215
216
217
218
    # 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
219
        conjugated = mrs.check_FID(repair=True)
220
221
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
222
223
                warnings.warn('FID has been checked and conjugated.'
                              ' Please check!', UserWarning)
224
225
226
227
228

    if args.conjbasis is not None:
        if args.conjbasis:
            mrs.conj_Basis()
    else:
William Clarke's avatar
William Clarke committed
229
        conjugated = mrs.check_Basis(repair=True)
230
231
        if args.verbose:
            if conjugated == 1:
William Clarke's avatar
William Clarke committed
232
233
234
                warnings.warn('Basis has been checked and conjugated.'
                              ' Please check!', UserWarning)

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

239
    # Keep/Ignore metabolites
240
    mrs.keep(args.keep)
William Clarke's avatar
William Clarke committed
241
    mrs.ignore(args.ignore)
242

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

William Clarke's avatar
William Clarke committed
249
    ppmlim = args.ppmlim
250
    if ppmlim is not None:
William Clarke's avatar
William Clarke committed
251
        ppmlim = (ppmlim[0], ppmlim[1])
252
253

    # Parse metabolite groups
William Clarke's avatar
William Clarke committed
254
    metab_groups = misc.parse_metab_groups(mrs, args.metab_groups)
255
256

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

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

Saad Jbabdi's avatar
Saad Jbabdi committed
280
281
282
283
284
    if args.verbose:
        print(mrs)
        print('Fitting args:')
        print(Fitargs)

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

287
288
289
    # Quantification
    # Echo time
    if args.TE is not None:
William Clarke's avatar
William Clarke committed
290
        echotime = args.TE * 1E-3
William Clarke's avatar
William Clarke committed
291
    elif 'EchoTime' in FID.hdr_ext:
292
        echotime = FID.hdr_ext['EchoTime']
293
294
    else:
        echotime = None
295
296
297
298
299
300
301
    # Repetition time
    if args.TR is not None:
        repetition_time = args.TR
    elif 'RepetitionTime' in FID.hdr_ext:
        repetition_time = FID.hdr_ext['RepetitionTime']
    else:
        repetition_time = None
William Clarke's avatar
William Clarke committed
302

303
    # Internal and Water quantification if requested
304
    if (mrs.H2O is None) or (echotime is None) or (repetition_time is None):
305
        if echotime is None:
William Clarke's avatar
William Clarke committed
306
307
308
            warnings.warn('H2O file provided but could not determine TE:'
                          ' no absolute quantification will be performed.',
                          UserWarning)
309
310
311
312
313
        if repetition_time is None:
            warnings.warn('H2O file provided but could not determine TR:'
                          ' no absolute quantification will be performed.',
                          UserWarning)
        res.calculateConcScaling(mrs, internal_reference=args.internal_ref, verbose=args.verbose)
314
    else:
315
316
317
318
        # Form quantification information
        q_info = quantify.QuantificationInfo(echotime,
                                             repetition_time,
                                             mrs.names,
319
320
321
322
                                             mrs.centralFrequency / 1E6,
                                             water_ref_metab=args.wref_metabolite,
                                             water_ref_metab_protons=args.ref_protons,
                                             water_ref_metab_limits=args.ref_int_limits)
323
324
325
326
327
328

        if args.tissue_frac:
            q_info.set_fractions(args.tissue_frac)
        if args.h2o_scale:
            q_info.add_corr = args.h2o_scale

329
        res.calculateConcScaling(mrs,
330
331
332
                                 quant_info=q_info,
                                 internal_reference=args.internal_ref,
                                 verbose=args.verbose)
333

334
335
336
    # Combine metabolites.
    if args.combine is not None:
        res.combine(args.combine)
Saad Jbabdi's avatar
Saad Jbabdi committed
337
338
339
340
    stop = time.time()

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

William Clarke's avatar
William Clarke committed
347
348
349
350
351
352
353
354
    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')
355
    if args.algo == 'MH':
William Clarke's avatar
William Clarke committed
356
357
358
359
360
        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
361

Saad Jbabdi's avatar
Saad Jbabdi committed
362
    # Save image of MRS voxel
363
    location_fig = None
364
365
366
367
368
369
    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
370
371

    # Save quick summary figure
William Clarke's avatar
William Clarke committed
372
373
374
375
    report.fitting_summary_fig(mrs, res,
                               filename=os.path.join(args.output,
                                                     'fit_summary.png'))

Saad Jbabdi's avatar
Saad Jbabdi committed
376
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
377
    if args.report:
William Clarke's avatar
William Clarke committed
378
379
380
381
382
383
384
385
386
387
        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
388
389
390
    if args.verbose:
        print('\n\n\nDone.')

William Clarke's avatar
William Clarke committed
391
392
393

def str_or_int_arg(x):
    try:
394
        return int(x)
William Clarke's avatar
William Clarke committed
395
    except ValueError:
396
        return x
397

William Clarke's avatar
William Clarke committed
398

399
def tissue_frac_arg(x):
Saad Jbabdi's avatar
Saad Jbabdi committed
400
    import json
401
402
    try:
        with open(x) as jsonFile:
William Clarke's avatar
William Clarke committed
403
            jsonString = jsonFile.read()
404
        return json.loads(jsonString)
William Clarke's avatar
William Clarke committed
405
    except IOError:
406
407
        return float(x)

William Clarke's avatar
William Clarke committed
408

409
410
411
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
412
413
        if isinstance(values[0], dict):
            setattr(namespace, self.dest, values[0])
414
        else:
William Clarke's avatar
William Clarke committed
415
416
417
418
            setattr(namespace, self.dest,
                    {'WM': values[0], 'GM': values[1], 'CSF': values[2]})


Saad Jbabdi's avatar
Saad Jbabdi committed
419
420
if __name__ == '__main__':
    main()