The brainloc package serves to describe a location on a brain surface mesh that is in MNI305 space, like the surfaces of the fsaverage standard brain templates used by FreeSurfer. Typical use cases of this package in the field of computational neuroimaging include:
Statistical results are often presented as significant clusters of activation on a brain surface mesh (sets of connected vertices). Typical use cases for working with clusters are:
This vignette describes how to perform these tasks with the brainloc package.
You will need a directory with the brain templates you want to use, typically these are the fsaverage subjects that come with the FreeSurfer neuroimaging software suite. There are 2 ways to get them:
When FreeSurfer is installed and configured properly on your system, there is nothing you need to do. In that case, the environment variable SUBJECTS_DIR
points at the subjects
directory under the installation directory (FREESURFER_HOME
), and it will be used to find the templates by the get_subjects_dir
function. You can also set SUBJECTS_DIR
to your custom template directory, of course, and that will be used instead.
If you cannot or do not want to set environment variables, you can set the option brainloc.subjects_dir
in R to any directory that contains the expected directory structure like this: options("brainloc.subjects_dir"="~/my_fs_home/subjects/")
. This would mean that ~/my_fs_home/subjects/
exists on your system and holds the directories of the template subjects. We assume that you have the fsaverage
subject in there. Setting this option takes precedence over any environment variables.
This requires that you have an internet connection and accept the FreeSurfer software license. See the code in the first example of this vignette below, which downloads the templates. The templates will be downloaded only if they cannot be found, and downloaded ones will be found on the next call to the function, of course.
Let's say we have a coordinate in MNI305 surface space and want to find the closest vertex. We first load brainloc
and create a brainparc
object, which internally is a named list holding the surfaces (the left and right hemisphere meshes) and the surface parcellations (according to some brain atlas) for these meshes. The parcellations assign to each vertex of a mesh a single region.
One can construct a brainparc manually from arbitrary loaded meshes and parcellations using the brainparc
function, but here we use the function brainparc_fs
, which is more convenient if you want to read it from files in a SUBJECTS_DIR
of recon-all
output data. We create a brainparc for the fsaverage
template:
library("brainloc"); # Searches for FreeSurfer SUBJECTS_DIR on machine and downloads if needed. sjd = get_subjects_dir(allow_download = TRUE, accept_freesurfer_license = TRUE); bp = brainparc_fs(sjd, "fsaverage", surface = "white", atlas = "aparc");
By default, many functions in the brainloc
package print information to STDOUT, which is handy for creating command line clients from the library. We do not want that here, so we disable it globally here:
options("brainloc.silent"=TRUE);
Note: One can also decide for each function that prints something, they all have a silent
parameter.
Let's say we have some query coordinates in surface space and want to find the closest vertex in the brainparc
for each of the coordinates. Here is how we can do that:
query_coords = matrix(seq(9)+0.1, ncol = 3, byrow = TRUE); query_coords ccv = coord_closest_vertex(query_coords, get_surface(bp)); ccv
This assumes that the query coordinates are in the surfaces space and uses Euclidean distance. It works for all subjects and meshes.
This assumes that you query vertex (or coordinate) is in MNI305 space. This is true for all FreeSurfer standard template subjects, like fsaverage
, fsaverage6
, etc. Our brainparc
is for fsaverage
, so we can use the coordinates of its vertices to query.
query_vertices = c(10L, 145029L); query_hemis = c("lh", "rh"); # vertex 10 on left hemi and 145029L on right hemi. query_coords = get_surface_coords(bp, query_vertices, query_hemis); query_coords coord_info = coord_MNI305_info(query_coords, surface = get_surface(bp)); coord_info
Note: If you have the regfusionr
package installed, have a look at the method
parameter of the coord_MNI305_info
function: you can use the more accurate regfusionr
method to transform from MNI305 to MNI152 then. Open the docs with ?coord_MNI305_info
to learn more.
We can easily find the regions a vertex is assigned to in all atlases of a brain parcellation:
query_vertices = c(10L, 145029L); query_hemis = c("lh", "rh"); # vertex 10 on left hemi and 145029L on right hemi. regions = vertex_regions(bp, query_vertices, query_hemis); regions
One can also compute the distance of a vertex to all atlas regions. Several definitions for the distance of a vertex to a region are available, here we go with the defaults:
query_vertices = c(10L, 145029L); query_hemis = c("lh", "rh"); # vertex 10 on left hemi and 145029L on right hemi. region_dists = vertex_closest_regions(bp, query_vertices, query_hemis, num_regions_to_report=5L); region_dists
One can change the distance definition with the parameters distance
and linkage
of the vertex_closest_regions
function. Open the docs with ?vertex_closest_regions
to learn more.
A cluster is a set of vertices on a surface, and the statistical values that are assigned to the cluster's vertices. Clusters are used in vertex-wise analyses to define activated brain regions (where the exact meaning of activated depends on the context). A clusters has a freeform label or name, but typically only consecutive numbers are used.
The data structure used to store clusters in brainloc
is called clusterinfo
. The following two concepts are relevant to understanding how a clusterinfo
instance is typically constructed:
cluster overlay
is an integer vector that assigns to each surface vertex a cluster number (zero is assigned to background/non-cluster vertices).statistical map
contains arbitrary values, one for each surface vertex (a.k.a. per-vertex data). A thresholded map only contains values for the cluster vertices, the values for all non-cluster vertices are set to an arbitrary value (typically 0.0).A clusterinfo
object can be constructed for a subject in two ways:
cluster overlay
and the statistical map
(thresholded or unthresholded), use the clusterinfo
function.statistical map
, use the clusterinfo_from_thresholded_overlay
function instead.In the following example, we construct a clusterinfo
instance from a cluster overlay and the statistical map:
lh_tmap_file = system.file("extdata", "lh.tmap.mgh", package = "brainloc", mustWork = TRUE); rh_tmap_file = system.file("extdata", "rh.tmap.mgh", package = "brainloc", mustWork = TRUE); lh_overlay_file = system.file("extdata", "lh.cluster.overlayID.mgh", package = "brainloc", mustWork = TRUE); rh_overlay_file = system.file("extdata", "rh.cluster.overlayID.mgh", package = "brainloc", mustWork = TRUE); clinfo = clusterinfo(lh_overlay_file, rh_overlay_file, lh_tmap_file, rh_tmap_file, template_subject = "fsaverage", subjects_dir = sjd); is.clusterinfo(clinfo);
One can extract the clusters from the `clusterinfo`` instance like this:
clusters = get_clusters(clinfo);
This returns a named list, where the keys are cluster names and the values are integer vectors containing the vertices that belong the the respective cluster. See ?get_clusters
for more options, like returning only the clusters for a specific hemisphere.
The clusterinfo
object can now be used to compute various information on the clusters. We will start by computing information on the location of the clusters:
cl_details = cluster_location_details(clinfo); cl_details;
Note that the clusterinfo
object contains a brain parcellation (a brainparc
object, see at the top). We can also compute the overlap of the cluster with its atlas regions:
overlap = cluster_region_overlap(clinfo); overlap;
Sometimes one wants to find all peaks of a cluster. A peak exists at vertices which have the most extreme value in their neighborhood. For a brain mesh, the neighborhood of a vertex is defined as all vertices which are reachable by one hop along the mesh edges.
peaks = cluster_peaks(clinfo); peaks
This concludes the package vignette, thanks for reading.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.