featanalysis.py 17.6 KB
Newer Older
1
2
#!/usr/bin/env python
#
3
# featanalysis.py - Utility functions for loading/querying the contents of
4
5
6
7
8
# a FEAT analysis directory.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module provides a few utility functions for loading/querying the
9
10
11
12
13
14
15
contents of a FEAT analysis directory. They are primarily for use by the
:class:`.FEATImage` class, but are available for other uses if needed. The
following functions are provided:

.. autosummary::
   :nosignatures:

16
   isFEATImage
17
   isFEATDir
18
   hasStats
19
   hasMelodicDir
20
21
   getAnalysisDir
   getTopLevelAnalysisDir
22
   isFirstLevelAnalysis
23
24
   loadDesign
   loadContrasts
Paul McCarthy's avatar
Paul McCarthy committed
25
   loadFsf
26
27
28
29
30
31
32
33
34
35
36
37
   loadSettings
   getThresholds
   loadClusterResults


The following functions return the names of various files of interest:

.. autosummary::
   :nosignatures:

   getDataFile
   getResidualFile
38
   getMelodicFile
39
40
41
42
   getPEFile
   getCOPEFile
   getZStatFile
   getClusterMaskFile
43
44
45
"""


46
47
48
49
import                   collections
import                   logging
import os.path        as op
import numpy          as np
50

51
52
53
54
import fsl.utils.path       as fslpath
import fsl.transform.affine as affine
from . import image         as fslimage
from . import                  featdesign
55

56

57
58
59
log = logging.getLogger(__name__)


60
61
62
def isFEATImage(path):
    """Returns ``True`` if the given path looks like it is the input data to
    a FEAT analysis, ``False`` otherwise.
63
    """
64
65
66

    try:
        path = fslimage.addExt(path, mustExist=True)
Paul McCarthy's avatar
Paul McCarthy committed
67
    except fslimage.PathError:
68
69
        return False

70
    dirname  = op.dirname( path)
71
    filename = op.basename(path)
72

73
    return filename.startswith('filtered_func_data') and isFEATDir(dirname)
74
75


76
77
78
79
def isFEATDir(path):
    """Returns ``True`` if the given path looks like a FEAT directory, or
    looks like the input data for a FEAT analysis, ``False`` otherwise.  A
    FEAT directory:
80

81
82
83
84
85
     - Must be named ``*.feat``.
     - Must contain a file called ``filtered_func_data.nii.gz``.
     - Must contain a file called ``design.fsf``.
     - Must contain a file called ``design.mat``.
     - Must contain a file called ``design.con``.
86
87
88

    :arg path: A file / directory path.
    """
89

90
    path = op.abspath(path)
91

92
93
    if op.isdir(path): dirname = path
    else:              dirname = op.dirname(path)
94

95
96
    if not dirname.endswith('.feat'):
        return False
97

98
    try:
99
        fslimage.addExt(op.join(dirname, 'filtered_func_data'), mustExist=True)
100
    except fslimage.PathError:
101
        return False
102

103
104
105
    if not op.exists(op.join(dirname, 'design.fsf')): return False
    if not op.exists(op.join(dirname, 'design.mat')): return False
    if not op.exists(op.join(dirname, 'design.con')): return False
106

107
    return True
108
109


110
111
112
113
114
def hasStats(featdir):
    """Returns ``True`` if it looks like statistics have been calculated
    for the given FEAT analysis, ``False`` otherwise.
    """

115
116
117
    try:
        getZStatFile(featdir, 0)
        return True
Paul McCarthy's avatar
Paul McCarthy committed
118
    except fslimage.PathError:
119
        return False
120

121

122
123
124
125
126
127
128
def hasMelodicDir(featdir):
    """Returns ``True`` if the data for the given FEAT directory has had
    MELODIC run on it, ``False`` otherwise.
    """
    return op.exists(getMelodicFile(featdir))


129
130
131
132
def getAnalysisDir(path):
    """If the given path is contained within a FEAT directory, the path
    to that FEAT directory is returned. Otherwise, ``None`` is returned.
    """
133
    featdir = fslpath.deepest(path, ['.feat'])
134

135
136
    if featdir is not None and isFEATDir(featdir):
        return featdir
137

138
    return None
139
140
141


def getTopLevelAnalysisDir(path):
142
143
144
    """If the given path is contained within a hierarchy of FEAT or MELODIC
    directories, the path to the highest-level (i.e. the shallowest in the
    file system) directory is returned. Otherwise, ``None`` is returned.
145
    """
146
    return fslpath.shallowest(path, ['.ica', '.gica', '.feat', '.gfeat'])
147
148


149
150
151
152
def getReportFile(featdir):
    """Returns the path to the FEAT report index file, or ``None`` if there
    is no report.
    """
153

154
155
156
157
158
    report = op.join(featdir, 'report.html')
    if op.exists(report): return report
    else:                 return None


159
def loadContrasts(featdir):
160
    """Loads the contrasts from a FEAT directory. Returns a tuple containing:
161

162
      - A list of names, one for each contrast.
163

164
      - A list of contrast vectors (each of which is a list itself).
165
166

    :arg featdir: A FEAT directory.
167
168
169
170
171
172
    """

    matrix       = None
    numContrasts = 0
    names        = {}
    designcon    = op.join(featdir, 'design.con')
173
174

    log.debug('Loading FEAT contrasts from {}'.format(designcon))
175

176
177
178
179
180
181
182
183
184
    with open(designcon, 'rt') as f:

        while True:
            line = f.readline().strip()

            if line.startswith('/ContrastName'):
                tkns       = line.split(None, 1)
                num        = [c for c in tkns[0] if c.isdigit()]
                num        = int(''.join(num))
185

186
                # The /ContrastName field may not
187
188
189
190
                # actually have a name specified
                if len(tkns) > 1:
                    name       = tkns[1].strip()
                    names[num] = name
191
192
193
194
195
196
197

            elif line.startswith('/NumContrasts'):
                numContrasts = int(line.split()[1])

            elif line == '/Matrix':
                break

198
        matrix = np.loadtxt(f, ndmin=2)
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218

    if matrix       is None             or \
       numContrasts != matrix.shape[0]:
        raise RuntimeError('{} does not appear to be a '
                           'valid design.con file'.format(designcon))

    # Fill in any missing contrast names
    if len(names) != numContrasts:
        for i in range(numContrasts):
            if i + 1 not in names:
                names[i + 1] = str(i + 1)

    names     = [names[c + 1] for c in range(numContrasts)]
    contrasts = []

    for row in matrix:
        contrasts.append(list(row))

    return names, contrasts

Paul McCarthy's avatar
Paul McCarthy committed
219

220
def loadFsf(designfsf):
Paul McCarthy's avatar
Paul McCarthy committed
221
222
    """Loads the analysis settings from a text file (.fsf) used to configure
    FEAT.
223

224
    Returns a dict containing the settings specified in the file
225

226
    :arg designfsf: A .fsf file.
227
    """
Paul McCarthy's avatar
Paul McCarthy committed
228

229
    settings  = collections.OrderedDict()
Paul McCarthy's avatar
Paul McCarthy committed
230

231
232
    log.debug('Loading FEAT settings from {}'.format(designfsf))

233
234
235
236
237
238
239
240
241
242
243
    with open(designfsf, 'rt') as f:

        for line in f.readlines():
            line = line.strip()

            if not line.startswith('set '):
                continue

            tkns = line.split(None, 2)

            key = tkns[1].strip()
244
            val = tkns[2].strip(' \'"')
245
246
247

            if key.startswith('fmri(') and key.endswith(')'):
                key = key[5:-1]
248

249
            settings[key] = val
250

251
252
    return settings

Paul McCarthy's avatar
Paul McCarthy committed
253

254
255
256
257
258
259
260
261
def loadSettings(featdir):
    """Loads the analysis settings from a FEAT directory.

    Returns a dict containing the settings specified in the ``design.fsf``
    file within the directory

    :arg featdir: A FEAT directory.
    """
Paul McCarthy's avatar
Paul McCarthy committed
262

263
    designfsf = op.join(featdir, 'design.fsf')
Paul McCarthy's avatar
Paul McCarthy committed
264

265
266
    return loadFsf(designfsf)

267

268
269
270
271
272
273
274
def loadDesign(featdir, settings):
    """Loads the design matrix from a FEAT directory.

    :arg featdir:  A FEAT directory.

    :arg settings: Dictionary containing FEAT settings (see
                   :func:`loadSettings`).
275

276
277
278
279
280
281
    :returns:      a :class:`.FEATFSFDesign` instance which represents the
                   design matrix.
    """
    return featdesign.FEATFSFDesign(featdir, settings)


282
def getThresholds(settings):
283
284
285
286
287
288
    """Given a FEAT settings dictionary, returns a dictionary of
    ``{stat : threshold}`` mappings, containing the thresholds used
    in the FEAT statistical analysis.

    The following keys will be present. Threshold values will be ``None``
    if the respective statistical thresholding was not carried out:
289

290
291
292
293
294
      - ``p``: P-value thresholding
      - ``z``: Z-statistic thresholding

    :arg settings: A FEAT settings dictionary (see :func:`loadSettings`).
    """
295
296
297
298
299
300
301
    p = settings.get('prob_thresh', None)
    z = settings.get('z_thresh',    None)

    if p is not None: p = float(p)
    if z is not None: z = float(z)

    return {'p' : p, 'z' : z}
302
303


304
305
306
307
308
309
310
311
312
def isFirstLevelAnalysis(settings):
    """Returns ``True`` if the FEAT analysis described by ``settings``
    is a first level analysis, ``False`` otherwise.

    :arg settings: A FEAT settings dictionary (see :func:`loadSettings`).
    """
    return settings['level'] == '1'


313
314
def loadClusterResults(featdir, settings, contrast):
    """If cluster thresholding was used in the FEAT analysis, this function
315
316
    will load and return the cluster results for the specified (0-indexed)
    contrast number.
317

318
319
    If there are no cluster results for the given contrast, ``None`` is
    returned.
320
321

    An error will be raised if the cluster file cannot be parsed.
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351

    :arg featdir:  A FEAT directory.
    :arg settings: A FEAT settings dictionary.
    :arg contrast: 0-indexed contrast number.

    :returns:      A list of ``Cluster`` instances, each of which contains
                   information about one cluster. A ``Cluster`` object has the
                   following attributes:

                     ============ =========================================
                     ``index``    Cluster index.
                     ``nvoxels``  Number of voxels in cluster.
                     ``p``        Cluster p value.
                     ``logp``     :math:`-log_{10}` of the cluster P value.
                     ``zmax``     Maximum Z value in cluster.
                     ``zmaxx``    X voxel coordinate of maximum Z value.
                     ``zmaxy``    Y voxel coordinate of maximum Z value.
                     ``zmaxz``    Z voxel coordinate of maximum Z value.
                     ``zcogx``    X voxel coordinate of cluster centre of
                                  gravity.
                     ``zcogy``    Y voxel coordinate of cluster centre of
                                  gravity.
                     ``zcogz``    Z voxel coordinate of cluster centre of
                                  gravity.
                     ``copemax``  Maximum COPE value in cluster.
                     ``copemaxx`` X voxel coordinate of maximum COPE value.
                     ``copemaxy`` Y voxel coordinate of maximum COPE value.
                     ``copemaxz`` Z voxel coordinate of maximum COPE value.
                     ``copemean`` Mean COPE of all voxels in the cluster.
                     ============ =========================================
352
353
    """

354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    # Cluster files are named like
    # 'cluster_zstatX.txt', where
    # X is the COPE number. And
    # the ZMax/COG etc coordinates
    # are usually in voxel coordinates
    coordXform  = np.eye(4)
    clusterFile = op.join(
        featdir, 'cluster_zstat{}.txt'.format(contrast + 1))

    if not op.exists(clusterFile):

        # If the analysis was performed in standard
        # space (e.g. a higher level group analysis),
        # the cluster file will instead be called
        # 'cluster_zstatX_std.txt', so we'd better
        # check for that too.
        clusterFile = op.join(
            featdir, 'cluster_zstat{}_std.txt'.format(contrast + 1))

373
374
375
        if not op.exists(clusterFile):
            return None

376
        # In higher level analysis run in some standard
377
378
379
380
381
        # space, the cluster coordinates are in standard
        # space. We transform them to voxel coordinates.
        # later on.
        coordXform = fslimage.Image(
            getDataFile(featdir),
382
            loadData=False).worldToVoxMat
383
384
385
386
387
388
389
390
391
392
393

    log.debug('Loading cluster results for contrast {} from {}'.format(
        contrast, clusterFile))

    # The cluster.txt file is converted
    # into a list of Cluster objects,
    # each of which encapsulates
    # information about one cluster.
    class Cluster(object):
        def __init__(self, **kwargs):
            for name, val in kwargs.items():
394

395
                attrName = colmap[name]
396
                if val is not None:
397
                    val = float(val)
398

399
400
                setattr(self, attrName, val)

401
402
403
404
405
406
            # if cluster thresholding was not used,
            # the cluster table will not contain
            # P valuse.
            if not hasattr(self, 'p'):    self.p    = 1.0
            if not hasattr(self, 'logp'): self.logp = 0.0

407
    # This dict provides a mapping between
408
409
    # Cluster object attribute names, and
    # the corresponding column name in the
410
    # cluster.txt file.
411
    colmap = {
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
        'Cluster Index'    : 'index',
        'Voxels'           : 'nvoxels',
        'P'                : 'p',
        '-log10(P)'        : 'logp',
        'Z-MAX'            : 'zmax',
        'Z-MAX X (vox)'    : 'zmaxx',
        'Z-MAX Y (vox)'    : 'zmaxy',
        'Z-MAX Z (vox)'    : 'zmaxz',
        'Z-COG X (vox)'    : 'zcogx',
        'Z-COG Y (vox)'    : 'zcogy',
        'Z-COG Z (vox)'    : 'zcogz',
        'Z-MAX X (mm)'     : 'zmaxx',
        'Z-MAX Y (mm)'     : 'zmaxy',
        'Z-MAX Z (mm)'     : 'zmaxz',
        'Z-COG X (mm)'     : 'zcogx',
        'Z-COG Y (mm)'     : 'zcogy',
        'Z-COG Z (mm)'     : 'zcogz',
        'COPE-MAX'         : 'copemax',
        'COPE-MAX X (vox)' : 'copemaxx',
        'COPE-MAX Y (vox)' : 'copemaxy',
        'COPE-MAX Z (vox)' : 'copemaxz',
        'COPE-MAX X (mm)'  : 'copemaxx',
        'COPE-MAX Y (mm)'  : 'copemaxy',
        'COPE-MAX Z (mm)'  : 'copemaxz',
        'COPE-MEAN'        : 'copemean'}
437
438
439
440
441
442
443
444
445
446
447
448

    # An error will be raised if the
    # cluster file does not exist (e.g.
    # if the specified contrast index
    # is invalid)
    with open(clusterFile, 'rt') as f:

        # Get every line in the file,
        # removing leading/trailing
        # whitespace, and discarding
        # empty lines
        lines = f.readlines()
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482

    lines = [line.strip() for line in lines]
    lines = [line         for line in lines if line != '']

    # the first line should contain column
    # names, and each other line should
    # contain the data for one cluster
    colNames     = lines[0]
    clusterLines = lines[1:]

    # each line should be tab-separated
    colNames     =  colNames.split('\t')
    clusterLines = [cl      .split('\t') for cl in clusterLines]

    # Turn each cluster line into a
    # Cluster instance. An error will
    # be raised if the columm names
    # are unrecognised (i.e. not in
    # the colmap above), or if the
    # file is poorly formed.
    clusters = [Cluster(**dict(zip(colNames, cl))) for cl in clusterLines]

    # Make sure all coordinates are in voxels -
    # for first level analyses, the coordXform
    # will be an identity transform (the coords
    # are already in voxels). But for higher
    # level, the coords are in mm, and need to
    # be transformed to voxels.
    for c in clusters:

        zmax    = [c.zmaxx,    c.zmaxy,    c.zmaxz]
        zcog    = [c.zcogx,    c.zcogy,    c.zcogz]
        copemax = [c.copemaxx, c.copemaxy, c.copemaxz]

483
484
485
        zmax    = affine.transform([zmax],    coordXform)[0]
        zcog    = affine.transform([zcog],    coordXform)[0]
        copemax = affine.transform([copemax], coordXform)[0]
486

487
488
489
        c.zmaxx,    c.zmaxy,    c.zmaxz    = zmax
        c.zcogx,    c.zcogy,    c.zcogz    = zcog
        c.copemaxx, c.copemaxy, c.copemaxz = copemax
490
491

    return clusters
492
493


494
def getDataFile(featdir):
495
    """Returns the name of the file in the FEAT directory which contains
496
    the model input data (typically called ``filtered_func_data.nii.gz``).
497

498
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
499

500
    :arg featdir: A FEAT directory.
501
    """
502
503
    datafile = op.join(featdir, 'filtered_func_data')
    return fslimage.addExt(datafile, mustExist=True)
504
505


506
def getMelodicFile(featdir):
507
    """Returns the name of the file in the FEAT results which contains the
508
    melodic components (if melodic ICA was performed as part of the FEAT
509
510
    analysis). This file can be loaded as a :class:`.MelodicImage`.

511
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
512
    """
513
514
    melfile = op.join(featdir, 'filtered_func_data.ica', 'melodic_IC')
    return fslimage.addExt(melfile, mustExist=True)
515
516


517
518
519
def getResidualFile(featdir):
    """Returns the name of the file in the FEAT results which contains
    the model fit residuals (typically called ``res4d.nii.gz``).
520

521
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
522

523
    :arg featdir: A FEAT directory.
524
    """
525
526
    resfile = op.join(featdir, 'stats', 'res4d')
    return fslimage.addExt(resfile, mustExist=True)
527

528

529
def getPEFile(featdir, ev):
530
    """Returns the path of the PE file for the specified EV.
531

532
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
533

534
535
536
    :arg featdir: A FEAT directory.
    :arg ev:      The EV number (0-indexed).
    """
537
538
    pefile = op.join(featdir, 'stats', 'pe{}'.format(ev + 1))
    return fslimage.addExt(pefile, mustExist=True)
539
540
541


def getCOPEFile(featdir, contrast):
542
543
    """Returns the path of the COPE file for the specified contrast.

544
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
545

546
    :arg featdir:  A FEAT directory.
547
    :arg contrast: The contrast number (0-indexed).
548
    """
549
550
    copefile = op.join(featdir, 'stats', 'cope{}'.format(contrast + 1))
    return fslimage.addExt(copefile, mustExist=True)
551
552


553
def getZStatFile(featdir, contrast):
554
555
    """Returns the path of the Z-statistic file for the specified contrast.

556
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
557

558
    :arg featdir:  A FEAT directory.
559
    :arg contrast: The contrast number (0-indexed).
560
    """
561
562
    zfile = op.join(featdir, 'stats', 'zstat{}'.format(contrast + 1))
    return fslimage.addExt(zfile, mustExist=True)
563
564
565


def getClusterMaskFile(featdir, contrast):
566
567
    """Returns the path of the cluster mask file for the specified contrast.

568
    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
569

570
    :arg featdir:  A FEAT directory.
571
    :arg contrast: The contrast number (0-indexed).
572
    """
573
574
    mfile = op.join(featdir, 'cluster_mask_zstat{}'.format(contrast + 1))
    return fslimage.addExt(mfile, mustExist=True)