05-CompareBrainRegionsWithCustomParcellation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare brain regions with a custom parcellation\n",
    "\n",
    "It is often useful to determine spatial overlap or correlations between maps from different brain parcellations, or between a region of interest and a brain map. This tutorial shows a simple case, where a custom map of the Destrieux atlas in MNI152 space is assessed against a group of probabilistic region maps from the Julich-Brain. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import siibra\n",
    "from nilearn import plotting\n",
    "import numpy as np\n",
    "import re\n",
    "atlas = siibra.atlases['human']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a custom Destrieux parcellation from region names and a labelled nifti volume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Destrieux atlas is ported in Freesurfer and it is provided in the package for standard surfaces. To obtain a decent volumetric form, a transformation from fsaverage to MNI152 was performed using the registration fusion technique by Wu et al. (https://pubmed.ncbi.nlm.nih.gov/29770530/). \n",
    "\n",
    "The volumetric map was compute and kindly provided by C. Paquola (Forschungszentrum Julich) and hosted on a public owncloud share (https://fz-juelich.sciebo.de/s/6KwoaksYyhU3br6). Note that the volumetric map is not perfect, it contains some discontinuities. \n",
    "\n",
    "For the numbering of parcels, the FreesurferColorLUT.txt which is shipped with [freesurfer](https://surfer.nmr.mgh.harvard.edu/) can be used. The list for Destrieux regions starts at Line 1234, with their label indices having an offset of 10100 (eg: 1038 in the volume is given as 11138 in the lookup table). For convenience, a copy of this file is also available on the above shared folder. \n",
    "\n",
    "The first step is to retrieve the labelled volume and region metadata files from the owncloud share. `siibra` has a convenient helper class for grabbing files from this type of repository:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "share = siibra.retrieval.OwncloudConnector(\n",
    "    server = \"https://fz-juelich.sciebo.de\",\n",
    "    share = \"6KwoaksYyhU3br6\"\n",
    ")\n",
    "labelvolume = share.get(\"aparc.a2009s.mni152.nii\")\n",
    "regionspecs = share.get(\"FreeSurferColorLUT.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We extract labelindex/regionname mappings from the text file, and use them to create a new `Parcellation` object for the Destrieux atlas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the regionname and labelindex from text to define regions\n",
    "lines = regionspecs.split(\"\\n\")\n",
    "regiondefs = [\n",
    "    {\n",
    "        'labelIndex': int(l.split()[0])-10100,\n",
    "        'name': l.split()[1]\n",
    "    }\n",
    "    for l in lines[1233:] \n",
    "    if re.match(f'\\d.*$',l)\n",
    "]\n",
    "\n",
    "# Build parcellation object\n",
    "destrieux = siibra.core.Parcellation(\n",
    "    identifier=\"destrieux\",\n",
    "    name=\"Volumetric Destrieux parcellation\",\n",
    "    regiondefs = regiondefs,\n",
    "    maps = [\n",
    "        siibra.volumes.LocalNiftiVolume(\n",
    "            name = \"Destrieux parcellation in MNI152\",\n",
    "            img = labelvolume,\n",
    "            space = atlas.get_space('mni152'))]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing overlaps of specific probabilistic maps with a parcellation map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to compute correlations of specific regions of Julich-Brain with another parcellations. \n",
    "\n",
    "`siibra`'s `ParcellationMap` class provides the method `assign` to compute overlaps with a custom map, and returns a list of assignments that is sorted by correlation strength. Each assignment is a 3-tuple, containing the assigned region object, a reference to the used region mask, and a dictionary of overlap scores, which includes amongst others Pearson's correlation coefficient. Note that the first call to `assign` is quite slow, since it will build all the required regional masks. Subsequent calls are fast, since `siibra` caches the mask definitions. \n",
    "\n",
    "We are interested to match an arbitrary subset of Julich-Brain regions individually, as well as their combination into a \"group map\". This can be achieved by creating a custom group region from a regular expression search across region names, using `siibra`'s `find_regions` method.\n",
    "\n",
    "We put the extraction of the group region and matching to anther parcellation into a function so we can reuse it for different choices of regions and parcellations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign(jubrain_regex, other_map, min_correlation = 0.3):\n",
    "    \"\"\"\n",
    "    Assign Julich-brain regions matched by the given regex \n",
    "    to another parcellation map.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    min_correlation: float\n",
    "        only matches with a Pearson correlation stronger \n",
    "        than this value are considered.\n",
    "    \n",
    "    Return\n",
    "    ------\n",
    "        Dictionary of assignments per Julich-brain region. \n",
    "    \"\"\"\n",
    "    import re \n",
    "    assignments = {}\n",
    "\n",
    "    # build the group region in this hemisphere.\n",
    "    # It will include the group as a parent region, \n",
    "    # as well as the individual child regions.\n",
    "    groupregion, = atlas.find_regions(\n",
    "        re.compile(jubrain_regex), \n",
    "        build_groups=True, \n",
    "        groupname = \"Group \"+jubrain_regex)\n",
    "\n",
    "    # compute assignments across the group\n",
    "    for region in groupregion: #\n",
    "        \n",
    "        mask = region.build_mask(\"mni152\",maptype='continuous')        \n",
    "        result = other_map.assign(\n",
    "            mask, msg=f\"Assigning {region.name:20.20}\")\n",
    "        significant = [\n",
    "            r for r in result \n",
    "            if abs(r[2]['correlation']) > min_correlation ]\n",
    "        if len(significant)>0:\n",
    "            assignments[region] = significant\n",
    "                \n",
    "    return assignments\n",
    "                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also define a function for plotting resulting assignments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_assignments(assignments):\n",
    "    \n",
    "    for query_region, matches in assignments.items():\n",
    "        for (matched_region, matched_map, scores) in matches:\n",
    "            view = plotting.plot_stat_map(\n",
    "                query_region.build_mask(\"mni152\", maptype='continuous'),\n",
    "                title=f\"{matched_region.name} / {query_region.name} \"\n",
    "                    f\"(ccoeff={scores['correlation']:.2f})\")\n",
    "            view.add_contours(matched_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Match Destrieux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "destrieux_map = destrieux.get_map(space=\"mni152\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assign_destrieux_l = assign('Fo[4-7].*left', destrieux_map)\n",
    "plot_assignments(assign_destrieux_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assign_destrieux_r = assign('Fo[4-7].*right', destrieux_map)\n",
    "plot_assignments(assign_destrieux_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Match Desikan/Kiliany\n",
    "\n",
    "Instead of the custom Destrieux parcels, we perform the same analysis with the predefined DK parcellation shipped with `siibra`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dk_map = atlas.get_map('mni152','desikan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assignments_dk_l = assign('Fo[4-7].*left', dk_map)\n",
    "plot_assignments(assignments_dk_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assignments_dk_r = assign('Fo[4-7].*right', dk_map)\n",
    "plot_assignments(assignments_dk_r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}