searchlight_surf.ipynb

{
  "metadata": {
    "name": "Surface-based searchlight on fMRI data"
  }, 
  "nbformat": 3, 
  "nbformat_minor": 0, 
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "heading", 
          "level": 1, 
          "metadata": {}, 
          "source": [
            "Surface-based searchlight on fMRI data"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "This example employs a surface-based searchlight as described in\n", 
            "*Oosterhof et al. (2011)* (with a minor difference that distances\nare currently computed using a Dijkstra distance metric rather than a geodesic\none). For more details, see the ", 
            "[Surfing](http://surfing.sourceforge.net)\nwebsite.\n\n", 
            "Surfaces used in this example are available in the tutorial dataset files;\neither the tutorial_data_surf_minimal or tutorial_data_surf_complete version.\nThe surfaces were reconstructed using FreeSurfer and\nsubsequently preprocessed with AFNI and SUMA using the\npymvpa2-prep-afni-surf wrapper script in PyMVPA's 'bin' directory, which\nresamples the surfaces to standard topologies (with different resolutions)\nusing MapIcosehedron, aligns surfaces to a reference functional volume, and\nmerges left and right hemispheres into single surface files. A more detailed\ndescription of the steps that this script takes is provided in the\ndocumentation on the ", 
            "[Surfing](http://surfing.sourceforge.net)\nwebsite.\n\n", 
            "If you use the surface-based searchlight code for a publication, please cite\nboth ", 
            "*PyMVPA (2009)* and ", 
            "*Oosterhof et al. (2011)*.\n\n", 
            "As always, we first have to import PyMVPA."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "from mvpa2.suite import *\nfrom mvpa2.clfs.svm import LinearCSVMC\nfrom mvpa2.base.hdf5 import h5save, h5load"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "As searchlight analyses are usually quite expensive in term of computational\nresources, we are going to enable some progress output to entertain us while\nwe are waiting."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "if __debug__:\n    from mvpa2.base import debug\n    debug.active += [\"SVS\", \"SLC\"]"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Define surface and volume data paths:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "rootpath = os.path.join(pymvpa_datadbroot,\n                        'tutorial_data', 'tutorial_data')\n\ndatapath = os.path.join(rootpath, 'haxby2001')\nsurfpath = os.path.join(rootpath, 'suma_surfaces')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Define functional data volume filename:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "epi_fn = os.path.join(datapath, 'sub001', 'BOLD', 'task001_run001', 'bold.nii.gz')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In this example we are concerned with the left hemisphere only.\n(Other possible values are 'r' for the right hemisphere and 'm' for merged\nhemispheres; the latter contains the nodes from the left and right\nhemispheres in a single file. Both the 'r' and 'm' options require\nthe tutorial_data_surf_complete tutorial data.)"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "hemi = 'l'"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Define the surfaces that enclose the grey matter, which are used to\ndelineate the grey matter. The pial surface is the 'outside' border of the\ngrey matter; the white surface is the 'inside' border.\n\n", 
            "The surfaces in this example were resampled using AFNI's MapIcosahedron\n(for more details, see the top of this script). ld refers to the number of\nlinear divisions of the twenty 'large' triangles of the original\nicosahedron; ld=x means there are 10*x**2+2 nodes (a.k.a. vertices)\nand 20*x**2 triangles (a.k.a. faces)."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "highres_ld = 128 # 64 or 128 is reasonable\n\npial_surf_fn = os.path.join(surfpath, \"ico%d_%sh.pial_al.asc\"\n                                     % (highres_ld, hemi))\nwhite_surf_fn = os.path.join(surfpath, \"ico%d_%sh.smoothwm_al.asc\"\n                                      % (highres_ld, hemi))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Define the surface on which the nodes are centers of the searchlight. This\nsurface should be an 'intermediate' surface, which is formed by the\nnode-wise average spatial coordinates of the inner (white) and outer (pial)\nsurfaces.\n\n", 
            "In this example a surface coarser (fewer nodes) than the grey matter-enclosing\nsurfaces is employed. This reduces the number of searchlights and therefore\nthe script's execution time. Of course one could also use a surface that has\nthe same number of nodes as the grey-matter enclosing surfaces; this is\nactually the default and used when souce_surf_fn (assigned below) is set\nto None.\n\n", 
            "It is required that highres_ld is an integer multiple of lowres_ld, so that\nall nodes in the low-res surface have a corresponding node (i.e., with the\nsame, or almost the same, spatial coordinate) on the high-res surface.\n\n", 
            "Choice of lowres_ld and highres_ld is somewhat arbitrary and a trade-off\nbetween spatial specificity and execution speed. For highres_ld a value of at\nleast 64 is be advisable as this ensures enough anatomical detail is available\nto select voxels in the grey matter accurately. Typical values for lowres_ld\nrange from 8 to 64.\n\n", 
            "Note that the data in tutorial_data_surf_minimal only contains\nall necessary surfaces for visualization for lowres_ld=16. For other values\nof lowres_ld (4, 8, 32, 64 and 128) the surfaces in\ntutorial_data_surf_complete are required."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "lowres_ld = 16 # 16, 32 or 64 is reasonable. 4 and 8 are really fast\n\nsource_surf_fn = os.path.join(surfpath, \"ico%d_%sh.intermediate_al.asc\"\n                                             % (lowres_ld, hemi))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Radius is specified as either an int (referring to a fixed number of voxels\nacross searchlights, with a variable radius in millimeters (or whatever unit\nis used in the files that define the surfaces), or a float (referring to the\nradius in millimeters, with a variable number of voxels).\n\n", 
            "Note that \"a fixed number of voxels\" in this context actually means an\napproximation, in that on average that number of voxels is selected but the\nactual number will vary slightly (typically in the range +/- 2 voxels)"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "radius = 100"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We're all set to go to create a query engine that performs 'voxel\nselection', that is determines, for each node, which voxels are near it\n(that is, in the corresponding searchlight disc).\n\n", 
            "As a reminder, the only essential values we have set so far are the\nfilenames of three surfaces (high-res inner and outer,\nand low-res source surface), the functional volume, and the searchlight\nradius.\n\n", 
            "Note that if the functional data was preprocessed and subsequently masked,\nvoxel selection should take into account this mask. To do so, the\ninstantiation of the query engine below takes an optional argument\n'volume_mask' (which can be a PyMVPA dataset, a numpy array, a Nifti\nvolume, or a string representing the file name of a Nifti volume). It is,\nhowever, recommended to ", 
            "*not* mask the functional data prior to voxel\nselection, because the voxel selection uses (implicitly) a mask based on the\ngrey-matter enclosing surfaces already, and this mask is assumed to be more\nprecise than typical volume-based masking implementations.\n\n", 
            "Also note that, as described above, the argument defining the low-res source\nsurface can be omitted, in which case it is computed as the node-wise\naverage of the white and pial surface.)"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "qe = disc_surface_queryengine(radius, epi_fn,\n                              white_surf_fn, pial_surf_fn,\n                              source_surf_fn)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Voxel selection is now completed; each node has been assigned a list of\nlinear voxel indices in the searchlight. These result are stored in\n'qe.voxsel' and can be saved with h5save for later re-use.\n\n", 
            "(Linear voxel indices mean that each voxel is indexed by a value between\n0 (inclusive) and N (exclusive), where N is the number of voxels in the\nvolume (N = NX * NY * NZ, where NX, NY and NZ are the number of voxels in\nthe three spatial dimensions). For certain analyses one may want to index\nvoxels by 'sub indices' (triples (i,j,k) with 0<=i<NX, 0<=j<=NY,\nand 0<=k<NZ) or spatial coordinates; conversions amongst\nlinear and sub indices and spatial coordinates is provided by\nfunctions in the  VolGeom (volume geometry) instance stored in\n'qe.voxsel.volgeom'.)\n\n", 
            "From now on we follow the example as in doc/examples/searchlight.py.\n\n", 
            "First, cross-validation is defined using a (SVM) classifier."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "clf = LinearCSVMC()\n\ncv = CrossValidation(clf, NFoldPartitioner(),\n                     errorfx=lambda p, t: np.mean(p == t),\n                     enable_ca=['stats'])"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Set the roi_ids, that is the node indices that serve as searchlight\ncenter. In this example it is set to None, meaning that all nodes are used\nas a searchlight center. It is also possible to restrict the nodes that serve\nas a searchlight center: setting roi_ids=np.arange(400,800) means that only\nnodes in the range from 400 (inclusive) to 800 (exclusive) are used as a\nsearchlight center, and the result would be a partial brain map."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "roi_ids = None"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Combining the query-engine and the cross-validation defines the\nsearchlight. The postproc-step averages the classification accuracies\nin each cross-validation fold to a single overall classification accuracy.\n\n", 
            "Because roi_ids is None is this example it could be omitted - it is only\nincluded for instructive purposes."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "sl = Searchlight(cv, queryengine=qe, postproc=mean_sample(), roi_ids=roi_ids)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In the next step the functional data is loaded. We can reduce\nmemory requirements significantly by considering which voxels to load:\nsince the searchlight analysis will only use data from voxels that\nwere selected (at least once) by the voxel selection step, a mask is\nderived from the voxel selection results and used when loading the\nfunctional data."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "mask = qe.voxsel.get_mask()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Load the functional data for subject 1 and the condition model 1 in the\ndataset (object viewing with 8 object categories). Note that we use the\nmask that came from the voxel selection."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "model = 1\nsubject = 1\nof = OpenFMRIDataset(datapath)\ndataset = of.get_model_bold_dataset(model, subject,\n                                    noinfolabel='rest', mask=mask)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Apply some typical preprocessing steps"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "poly_detrend(dataset, polyord=1, chunks_attr='chunks')\n\ndataset = dataset[np.array([l in ['rest', 'house', 'scrambledpix']\n                           for l in dataset.targets], dtype='bool')]\n\nzscore(dataset, chunks_attr='chunks', param_est=('targets', ['rest']),\n        dtype='float32')\n\ndataset = dataset[dataset.sa.targets != 'rest']"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Run the searchlight on the dataset."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "sl_dset = sl(dataset)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Searchlight results are now stored in sl_dset. As sl_dset is just like\nany other PyMVPA dataset, it can be stored with h5save for future use.\n\n", 
            "The remainder of this example provides a data file that\ncan be visualized using AFNI's SUMA. This is achieved by storing the dataset\nas an NIML (NeuroImaging Markup Language) dataset that can be viewed by\nAFNI's SUMA. sl_dset contains a feature attribute 'center_ids' that is\nautomagically used to define the node indices of the searchlight centers in\nthis NIML dataset.\n\n", 
            "Note that this conversion will not preserve all information in sl_dset but\nonly the samples and (feature, sample, dataset) attributes that behave\nlike arrays or strings or scalars. For example, in this example sl_dset has a\ndataset attribute 'mapper' which is not stored in the NIML dataset (and\na warning message is printed during the conversion, which can be ignored\nsavely). As mentioned above, using h5save will preserve this information\n(but its output cannot be viewed in SUMA).\n\n", 
            "Before saving the dataset, first the labels are set for each sample (in\nthis case, only one) so that they show up in SUMA."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "sl_dset.sa['labels'] = ['HOUSvsSCRM']"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Set the filename for output.\nSearchlight results are stored in the surface directory for easy\nvisualization. Finally print an informative message on how the\ngenerated data can be visualized using SUMA."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# save as NIML dataset\nfn = 'ico%d-%d_%sh_%dvx.niml.dset' % (lowres_ld, highres_ld, hemi, radius)\npath_fn = os.path.join(surfpath, fn)\nniml.write(path_fn, sl_dset)\n\n# save as GIFTI\nif externals.exists('nibabel'):\n    fn = 'ico%d-%d_%sh_%dvx.func.gii' % (lowres_ld, highres_ld, hemi, radius)\n    path_fn = os.path.join(surfpath, fn)\n    map2gifti(sl_dset, path_fn)\n\n\nprint (\"To view results in SUMA, cd to '%s', run 'suma -spec \"\n      \"%sh_ico%d_al.spec', press ctrl+s, \"\n       \"click on 'Load Dset', and select %s\" %\n       (surfpath, hemi, lowres_ld, fn))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }
      ], 
      "metadata": {}
    }
  ]
}