Multivariate Pattern Analysis in Python

# Full-scale fMRI data analysis using pattern classification

This script demonstrates a complete classification analysis as it could be found in actual research projects. It is a representative example of the functionality accessible from the command line interface – when it was originally release with PyMVPA 2.3. We start by creating a dataset as a collection of information from various sources. fMRI data is loaded from a 4D NIfTI image, while only voxel with non-zero value in a mask image are kept. Each volume is associated with some attribute values that are read from a text file. In addition, a number of mask additional mask images are included as feature attributes, in order to be able to conveniently group voxel/features into ROIs. Lastly, motion estimates are included as per-volume attributes. The resulting dataset is stored in a compressed HDF5 file.

    pymvpa2 mkds \
       --openfmri-modelbold "$dataroot" 1 1 '25mm' \
       --mask "$dataroot"/sub001/masks/25mm/brain.nii.gz \
       --add-sa-attr "$dataroot"/../attributes_literal.txt  \
       --add-vol-attr hoc "$dataroot"/sub001/masks/25mm/hoc.nii.gz \
       --add-vol-attr gm "$dataroot"/sub001/masks/25mm/gray.nii.gz \
       --add-vol-attr vt "$dataroot"/sub001/masks/25mm/vt.nii.gz \
       --add-fsl-mcpar bold_moest.txt \
       --hdf5-compression gzip \
       -o "$outdir"/bold.hdf5

The describe command generates a terse summary of the dataset.

    pymvpa2 describe -i "$outdir"/bold.hdf5

One of the additional feature attributed was down-sampled and aligned brain parcelation of the Harvard-Oxford cortical atlas. The `dump` command can be used to extract any dataset component and convert it into a variety of formats – including plain text – for further processing with standard UNIX console tools.

    echo "Number of ROIs in the Harvard-Oxford cortial atlas: "
    pymvpa2 dump --fa hoc -f txt -i "$outdir"/bold.hdf5 | sort | uniq | wc -l

The `preproc` command offers a few selected preprocessing procedures. Here we first regress out the motion parameter estimate time-courses that were included as sample attributes in the dataset. And finally, we perform spectral filtering using a Butterworth with selected pass and stop bands. As with every processing step, the result is stored as a dataset in an HDF5 file.

    pymvpa2 preproc --poly-detrend 0 \
                    --detrend-regrs bold_moest.txt_0 bold_moest.txt_1 bold_moest.txt_2 bold_moest.txt_3 bold_moest.txt_4 bold_moest.txt_5 \
                    --filter-passband 0.005 0.067 \
                    --filter-stopband 0.0025 0.1 \
                    --sampling-rate 0.4 \
                    --chunks chunks \
                    --hdf5-compression gzip \
                    -o "$outdir"/bold_mcf.hdf5 \
                    -i "$outdir"/bold.hdf5

For the subsequent classification analysis we are only interested in two out of all available experimental conditions, and only in voxels that are part of the “VT” mask. Sub-selection of dataset content is supported by indices, as well as a simple expressions.

    pymvpa2 select \
        --samples-by-attr targets eq face or targets eq house \
        --features-by-attr vt gt 0 \
        -i "$outdir"/bold_mcf.hdf5 \
        -o "$outdir"/faceshouses_inVT.hdf5

The easiest way to perform a classification analysis is the selection of a pre-crafted classifier instance from the “warehouse”. A variety of data partitioning schemes is available – we select “leave-one-out” – by default operating on the `chunks` sample attribute of the dataset. Selection of error functions is possible by name, or by providing a Python script with a custom implementation. Many commands and options of the command line interface can be fed and extended with custom Python scripts.

    pymvpa2 crossval \
        --learner "$clf" \
        --partitioner n-1 \
        --errorfx mean_match_accuracy \
        --avg-datafold-results \
        -i "$outdir"/faceshouses_inVT.hdf5 \
        -o "$outdir"/xval_faces_vs_houses_inVT.hdf5

After this first ROI-based classification analysis, we are now aiming for a very similar classification analysis that, in contrast, is done in a “searchlight” – a traveling ROI analysis throught the entire brain. Hence we create a new dataset, again with only face and house data samples, but this time including all voxels.

    pymvpa2 select \
        --samples-by-attr targets eq face or targets eq house \
        --hdf5-compression gzip \
        -i "$outdir"/bold_mcf.hdf5 \
        -o "$outdir"/faceshouses_brain.hdf5

The `searchlight` command can be used to compute arbitrary metric in this fashion, but has built-in support for cross-validated classification analyses (`--payload`). Spherical ROI with a radius of 4 voxels (`--neighbors`) will be generated, centered on gray-matter voxels only (`--roi-attr`). The computation will be parallelized with up to two concurrent processes. [Note: The `--scatter-rois` option is only present to speed up computation

> and can be removed in order to obtain a dense result map.]

    pymvpa2 --dbg-channel SLC searchlight \
        --payload cv \
        --neighbors 4 \
        --scatter-rois 5 \
        --roi-attr gm \
        --nproc 2 \
        --cv-learner "$clf" \
        --cv-partitioner oddeven:chunks \
        --cv-errorfx mean_match_accuracy \
        --cv-avg-datafold-results \
        --cv-permutations 2 \
        --hdf5-compression gzip \
        -i "$outdir"/faceshouses_brain.hdf5 \
        -o "$outdir"/sl_faces_vs_houses_brain.hdf5

Using the `dump` command, results can be stored in various formats, including NIfTI. Saving as NIfTI automatically takes care of projecting back results into the 3D voxel space.

    pymvpa2 dump -s \
        -f nifti \
        -i "$outdir"/sl_faces_vs_houses_brain.hdf5 \
        -o "$outdir"/sl_faces_vs_houses_brain_ACC.nii.gz
    pymvpa2 dump --fa null_prob \
        -f nifti \
        -i "$outdir"/sl_faces_vs_houses_brain.hdf5 \
        -o "$outdir"/sl_faces_vs_houses_brain_NP.nii.gz

An alternative to a searchlight with its often arbitrary ROI shapes and boundaries is an iterative ROI analysis – cycling through a number of ROIs that are defined by localizers or an atlas. Here we perform the cross-validated classification analysis shown above on all areas defined in the Harvard-Oxford cortical atlas and present in our data.

    hoc_rois=( $(pymvpa2 exec -i "$outdir"/bold.hdf5 -e 'print(" ".join(map(str, dss[0].fa["hoc"].unique)))') )
    echo "ROIs of the Harvard-Oxford cortial atlas present in the data: ${hoc_rois[*]}"

    for roi in ${hoc_rois[*]}; do
        echo -en " ROI $roi\t"

Select corresponding voxels.

    pymvpa2 select --features-by-attr hoc eq $roi \
        -i "$outdir"/faceshouses_brain.hdf5 \
        -o "$outdir"/roi_tmp.hdf5

Report number of voxels present in the given ROI.

    nfeatures=$(pymvpa2 exec -i "$outdir"/roi_tmp.hdf5 -e "print(dss[0].nfeatures)")
    resultds="${outdir}/xval_faces_vs_houses_inROI${roi}.hdf5"

    echo -en "$nfeatures voxels\taccuracy="

And run the cross validation, finally printing the overall accuracy as a result on the console.

        pymvpa2 crossval \
            --learner "$clf" \
            --partitioner n-1 \
            --errorfx mean_match_accuracy \
            --avg-datafold-results \
            -i "$outdir"/roi_tmp.hdf5 \
            -o $resultds | awk '/ACC%/{printf "%.2f%%\n", $2}'

        [ -z "${MVPA_TESTS_QUICK:-}" ] || break  # reserved for testing
    done

See also

The full source code of this example is included in the PyMVPA source distribution (`doc/examples/cmdline/fmri_analyses.sh`).