[pymvpa] Surface searchlight problems (alignment related)

Nick Oosterhof n.n.oosterhof at googlemail.com
Mon Feb 9 15:10:28 UTC 2015

Hello Mike,

On 08 Feb 2015, at 20:58, Mike E. Klein <michaeleklein at gmail.com> wrote:

> I am attempting to do a surface-informed MVPA analysis and am approaching the “screw it, I’ll stick with my volumetric results” decision point

I’m sorry to hear that. I wrote most surface-based functionality in PyMVPA, so let’s see and try to resolve this.

> 		• 3. pymvpa2-pre-afni-surf run to:
> 			• a. resamples the surfaces to standard topologies (with different resolutions) using MapIcosehedron 
> 			• b. aligns surfaces to a reference functional volume, 
> 			• c. merges left and right hemispheres into single surface files.
> 			• I did this with pymvpa2-prep-afni-surf -e mySub/example_functional_volume.nii -d surfer/mySub/surf/ -r mySub/new_suma_surfaces/
> 				• example_functional_volume.nii is taken from my 4D time series.

Please note that the mySub/example_functional_volume.nii file is the one that is used for coregistration with the anatomical (and with the surfaces). If the EPI image is of poor quality, then coregistration may be poor as well. Also, typically the “-e” volume should be the same as the one used as reference for motion correction.

As an alternative, if you have coregistered an anatomical image to an EPI image yourself and are confident that the coregistration is good, you can also specify the anatomical volume with the “-a” option (and leave out the “-e” option).

> I visually verified that the entire time series was in good alignment with the example. 
> 		• Step 3a vs. Step 3b is a point of a confusion for me. I think I am being thrown off by the phrase “standard topologies.” We are projecting the anatomically-extracted surface into the functional space so that the functional data does not need to be transformed/interpolated (unintentionally smoothed). If that is the case, then how are “standard topologies” involved? (When I hear standard, I think “standard space” like MNI152, etc.) Is the idea that the -shape- of the surface sheet isn’t altered in this “resampl[ing],” but instead just the specific locations (and numbering schemes) for the vertices/edges? This is the only explanation that makes sense to me. 

Indeed. Surfaces from different participants almost always have a different number of vertices, and also the topology (which defines which vertices make up triangles) varies. The resampling ensures that surfaces from different participants have the same number of nodes, the same topology, and also that corresponding nodes are more or less (within the limits of surface-based coregistration across participants) at the same spatial location. 

> 		• 4. Visually check the alignment of one of the new suma surfaces (e.g. mh_ico32_al.spec) with the example functional volume using AFNI/SUMA. I am doing this by running:
> 			• afni -niml &
> 			• suma -spec mh_ico32_al.spec -sv example_functional_volume.nii & 
> 				• (I’m showing this with both a T2 and a T1 image registered the example functional in the T2’s 3mm space, as it’s easier to see.) These images are shown in two of the screenshots.
> 			• pressing “t” in SUMA to link together AFNI and SUMA
> 			• This, to me, looks OK, but not great. There is clearly much correspondence between the surfaces and the volume, although the STG/HG area (on right) looks particularly bad to me. 

I agree, the alignment is decent but not awesome. Did you use a recent version of PyMVPA? If so, there should be a couple of QA images generated using AFNI’s @AddEdge script; these are the *_ec* and *_ee* files. Can you maybe post / DropBox those so that we can see how the alignment looks there (between EPI and anatomical).

> 			• Separately, I get a warning message about a force switch from the original view to talairach view. I’m not sure why this would be, as all of these volumes and surfaces should be in the functional space…

Could it be that example_functional_volume.nii is considered in talairach view?

> 		• 5. Run the PyMVPA surface searchlight script. 
> 			• I’m doing this as a sanity check: decoding sound vs. silence.
> 			• I’m using the merged (m) hemisphere option and I’m using what I think are essentially default settings:
> 				• highres_ld = 128 
> 				• pial_surf_fn = os.path.join(surfpath, "ico%d_%sh.pial_al.asc % (highres_ld, hemi))
> 				• white_surf_fn = os.path.join(surfpath, "ico%d_%sh.smoothwm_al.asc” % (highres_ld, hemi))
> 				• lowres_ld = 32 # 16, 32 or 64 is reasonable. 4 and 8 are really fast
> 				• source_surf_fn = os.path.join(surfpath, "ico%d_%sh.intermediate_al.asc % (lowres_ld, hemi))
> 				• radius = 123
> 				• qe = disc_surface_queryengine(radius, epi_fn, white_surf_fn, pial_surf_fn, source_surf_fn, nproc=4)
> 				• clf = LinearCSVMC()
> 				• cv = CrossValidation(clf, NFoldPartitioner(), errorfx=lambda p, t: np.mean(p == t), enable_ca=['stats’])

That all looks fine to me.

> 		• 6. I’m then checking the results with AFNI/SUMA. 
> 			• While I do get very high (almost 100% decoding) in certain regions for sound/silence, the peaks are not along the STG/HG as they should be. Instead they are considerably posterior and superior (the SMG and surrounding). This doesn’t pass the smell test, as is also at odds with my volumetric results, which show maximal decoding in and around A1.

Indeed, that does not seem what you want. This could be caused by bad coregistration. Another option worth trying (although it does not explain the lack of a peak in A1) is a smaller radius (fewer voxels), to restrict the size of the searchlight.

> 			• A stray thought, but if I display the suma surfaces over my anatomical in its native space, the surface-HG region fits over the native anatomical IPL region. 

You mean that coregistration seems poor between surfaces and the SurfVol_al2exp+orig volume? That’s weird, because that coregistration should be just as good as the original coregistration between the input (FreeSurfer’s anatomical file and surface file, e.g. the stuff in the SUMA directory).

> 		• And, of course, even if I had all of the issues ironed out, I’m at a loss for how to treat the DSET surface MVPA files in a group analysis. In volume space, I’ve run my information maps through SPM. 

That’s a whole other issue, and there is little literature about doing this. I wrote functionality to do this in two ways:
1) the SPM way, i.e. computing residuals from the group mean and use those to estimate the smoothness. An experimental implementation is in PyMVPA/mvpa2/support/afni/afni_surface_alphasim.py 
2) Threshold-Free Cluster Enhancement (or any other statistic such as cluster size) using GNU Octave / Matlab. An example is here: http://cosmomvpa.org/contents_demo.html#demo-surface-tfce. Note that this requires (at the moment) NIML .dset files stored in *ASCII* format, which is not the default when using niml.write.

> At present, I can’t even get these searchlight result files to load in FreeSurfer’s FreeView (even if I convert file formats to something like GII).

That’s a bit surprising. Did you use ConvertDset, or another program? Have you tried just writing them in ACSII format? If this remains an issue, it may be good if you could make an offending dataset available.

So in summary, my first hunch is that coregistration between the functional data and the surfaces is not optimal. Can you try to use an anatomical file as input for pymvpa2-pre-afni-surf with the “-a” option after you’ve verified that the anatomical file is properly co-registered with the EPI image?

Please let me know how it goes, I hope we can get your surface-based analysis to work.


More information about the Pkg-ExpPsy-PyMVPA mailing list