.. AUTO-GENERATED FILE -- DO NOT EDIT!

.. _example_cmdline_fmri_analyses:

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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.]

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

  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.

.. code-block:: sh

      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.

.. code-block:: sh

      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.

.. code-block:: sh

      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
  

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