HIBALL-winterschool-2023.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "999772f5",
   "metadata": {},
   "source": [
    "# Instantiation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e71a4c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import siibra\n",
    "from packaging.version import Version\n",
    "assert Version(siibra.__version__) >= Version('1.0a14')\n",
    "\n",
    "import matplotlib\n",
    "from nilearn import plotting\n",
    "\n",
    "# ignore the following lines at this point - we just touch some objects to trigger data loading \n",
    "# while we still take the introduction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b46e900c",
   "metadata": {},
   "source": [
    "For the purpose of this tutorial, we request some objects here already to retrieve data while your tutor still does the introductory remarks :-) Please just run the cell and ignore for now. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3af4ed33",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_brain = siibra.parcellations.get('julich 3.1')\n",
    "julich_pmaps = julich_brain.get_map('mni152', 'statistical')\n",
    "julich_pmaps.fetch(region='4a left')\n",
    "siibra.warm_cache()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbc1c43d",
   "metadata": {},
   "source": [
    "# Part I: Accessing reference spaces, parcellations and regions\n",
    "\n",
    "\n",
    "## Instance tables of key concepts\n",
    "\n",
    "`siibra` is structured around the key concepts atlas, reference space, parcellation and parcellation region. Each of these concepts has a specific type. `siibra` comes with preconfigured instances of these concepts, which can be easily accessed via *instance tables*. When you load siibra for the first time, it pulls this preconfiguration information from our online repository and caches it on your computer. Therefore, `siibra` will be a bit slower when you use it for the first time (or after a version upgrade).\n",
    "\n",
    "Here is an overview of the key concepts:\n",
    "\n",
    "| Concept | siibra type | instance table | description |\n",
    "| :-- | :-- | :-- | :-- |\n",
    "| Atlases | Atlas | `siibra.atlases` | A collection of related reference spaces and parcellations, typically per species |\n",
    "| Reference spaces | Space | `siibra.spaces` | 3D coordinate systems of the brain |\n",
    "| Parcellations | Parcellation | `siibra.parcellations` | Different brain parcellations schemes with their region hierarchies |\n",
    "| regions | Region | - | Structures defined within a parcellation, each representing a subtree of a parcellation's hierarchy| \n",
    "\n",
    "**Note:** These concepts are just semantic objects - they mostly give names and relationships to atlas concepts. We will deal with parcellation maps in the next notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd71f02",
   "metadata": {},
   "outputs": [],
   "source": [
    "import siibra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3affd86",
   "metadata": {},
   "source": [
    "## Select parcellations from their instance table\n",
    "\n",
    "The instance table for parcellations is `siibra.parcellations`. To get an overview, we can simply print its elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38937f61",
   "metadata": {},
   "outputs": [],
   "source": [
    "siibra.parcellations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc9a7167",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(siibra.parcellations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1be34e2",
   "metadata": {},
   "source": [
    "To actually access an element from an instance table, there are several options:\n",
    "\n",
    " - by tab-completion on an interactive Python shell. This is very convenient, since it allows you to browse while typing.\n",
    " - by fuzzy keyword matching via the get() function\n",
    " - By using the index operator `[]` with keywords or numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bbee7cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "siibra.parcellations.JULICH_BRAIN_CYTOARCHITECTONIC_ATLAS_V3_0_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "895779b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "siibra.parcellations.get('bundles')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa3fe22d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is equivalent to the above\n",
    "siibra.parcellations['bundles']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "753a3cfd",
   "metadata": {},
   "source": [
    "## Your turn! Browse predefined reference spaces\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "Access the MNI 152 ICBM 2009c nonl asymmetric space.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81cf418",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "81ee3f7c",
   "metadata": {},
   "source": [
    "## Properties of parcellation objects\n",
    "\n",
    "The parcellations contain some metadata, such as name, description and related publication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35beb962",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_brain = siibra.parcellations.get('julich')\n",
    "print(julich_brain.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8a2e00d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(julich_brain.description)\n",
    "\n",
    "for p in julich_brain.publications:\n",
    "    print(\"\\n\" + p['citation'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f339ba8",
   "metadata": {},
   "source": [
    "## Brain regions\n",
    "\n",
    "Providing all brain regions through instance tables would not be helpful, since there are thousands of them. Regions are organized into hierarchical trees attached to a parcellation.  A parcellation object in siibra is actually a special type of region, namely the top node of a region tree with some additional metadata. \n",
    "\n",
    "Other regions may be the root of a more fine-grained subtree, or leafs without children.\n",
    "\n",
    "We can print the subtree of each region object using the `Region.tree2str` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3de8d3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(julich_brain.tree2str())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaeaa8bf",
   "metadata": {},
   "source": [
    "To access individual brain regions, we typically pass over a unique part of their name to the parcellation (or any parent region):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "220ff890",
   "metadata": {},
   "outputs": [],
   "source": [
    "amygdala = julich_brain.get_region('amygdala')\n",
    "print(amygdala.tree2str())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6d7a255",
   "metadata": {},
   "source": [
    "# Part II - Fetching reference templates and parcellation maps\n",
    "\n",
    "## Volumetric maps\n",
    "\n",
    "The same parcellation can be mapped in different refernce spaces. So to request a parcellation map, we have to link a parcellation with a reference space. In siibra, we request the map from the parcellation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4386e001",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_brain = (\n",
    "    siibra.parcellations.JULICH_BRAIN_CYTOARCHITECTONIC_ATLAS_V3_0_3\n",
    ")\n",
    "julich_mpm = julich_brain.get_map(\n",
    "    space=siibra.spaces.MNI_152_ICBM_2009C_NONLINEAR_ASYMMETRIC\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82b02dc6",
   "metadata": {},
   "source": [
    "At this point, we just have a parcellation map object. The image volume has not yet been loaded. To do so, we use the `fetch()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be81255e",
   "metadata": {},
   "outputs": [],
   "source": [
    "mpm_img = julich_mpm.fetch()\n",
    "mpm_img"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a80cc6a4",
   "metadata": {},
   "source": [
    "As you see, the image data is provided as a Nifti1Image, which is very common to use in neuroscience. siibra represents most image data in this format. We can use common libraries to visualize the image - we recommend the excellent `nilearn.plotting` tools. Some parcellation maps in siibra provide their own colormap that can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19bd1265",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_args = {\"symmetric_cmap\": False, \"colorbar\": False}\n",
    "plotting.plot_roi(\n",
    "    mpm_img,\n",
    "    title=julich_brain.name,\n",
    "    cmap=julich_mpm.get_colormap(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a84f4fd7",
   "metadata": {},
   "source": [
    "We can also fetch the map of a single region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d74eaea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "motor_map = julich_mpm.fetch(region='4a left')\n",
    "plotting.plot_roi(motor_map, title='4a left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c2a7667",
   "metadata": {},
   "source": [
    "## Regional probability maps\n",
    "\n",
    "So far we looked at labelled parcellation maps. To access probability maps, we explicitly request the \"statistical\" maptype. Everything else works the same way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dda37c8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_pmaps = julich_brain.get_map(\n",
    "    siibra.spaces.MNI_152_ICBM_2009C_NONLINEAR_ASYMMETRIC,\n",
    "    maptype=siibra.MapType.STATISTICAL\n",
    ")\n",
    "motor_pmap = julich_pmaps.fetch(region='4a left')\n",
    "plotting.plot_stat_map(motor_pmap, title='4a left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eb3fd45",
   "metadata": {},
   "source": [
    "## Surface maps\n",
    "\n",
    "Julich-Brain is also provided as a surface map in fsaverage. fsaverage represents a different space. This time, instead of a NIfTI image, we obtain a mesh structure, which is a dictionary with vertices, faces and vertex labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daf44d70",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_surfmap = julich_brain.get_map(siibra.spaces.FREESURFER_FSAVERAGE)\n",
    "surf = julich_surfmap.fetch()\n",
    "surf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "263b1875",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.view_surf(\n",
    "    surf_mesh=[surf['verts'], surf['faces']],\n",
    "    surf_map=surf['labels'], \n",
    "    cmap=julich_surfmap.get_colormap(),\n",
    "    **plot_args\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2fe1a26",
   "metadata": {},
   "source": [
    "### Your turn!\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "Load and display the DiFuMo 64 map in MNI152 space.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08e24aa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b3119b3",
   "metadata": {},
   "source": [
    "## Reference templates\n",
    "\n",
    "Just as we can fetch maps, we can fetch the reference template of a space. Let's do this to fetch the BigBrain model!\n",
    "\n",
    "Per default, fetch() will download a low-resolution version. We will talk about higher resolutions later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3016b1b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_brain.shortname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e8dba58",
   "metadata": {},
   "outputs": [],
   "source": [
    "bigbrain_space = siibra.spaces.BIGBRAIN_MICROSCOPIC_TEMPLATE_HISTOLOGY\n",
    "bigbrain_template = bigbrain_space.get_template()\n",
    "plotting.plot_img(bigbrain_template.fetch(), bg_img=None, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ddc323a",
   "metadata": {},
   "source": [
    "# Part III - Assigning locations to brain structures\n",
    "\n",
    "## Specifying brain locations\n",
    "\n",
    "Locations in the brain are often specified by coordinates, or peaks/blobs in images. siibra has predefined location types for points, sets of points, and bounding boxes in reference spaces.\n",
    "\n",
    "siibra tries to ensure that location objects, just as maps and images, are clearly linked to a reference space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1dc2620",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a point by specifying coordinates. You can do this in the viewer!\n",
    "point = siibra.locations.Point(\"24.150mm, -18.000mm, 42.150mm\", space='mni152')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9708f0c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_img(mpm_img, cut_coords=tuple(point), cmap=julich_mpm.get_colormap())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6da5d084",
   "metadata": {},
   "outputs": [],
   "source": [
    "point = siibra.from_json(\"\"\"\n",
    "{\n",
    "  \"@id\": \"676457d9\",\n",
    "  \"@type\": \"https://openminds.ebrains.eu/sands/CoordinatePoint\",\n",
    "  \"coordinateSpace\": {\n",
    "    \"@id\": \"minds/core/referencespace/v1.0.0/a1655b99-82f1-420f-a3c2-fe80fd4c8588\"\n",
    "  },\n",
    "  \"coordinates\": [\n",
    "    {\n",
    "      \"@id\": \"9837b595\",\n",
    "      \"@type\": \"https://openminds.ebrains.eu/core/QuantitativeValue\",\n",
    "      \"value\": 13.0557,\n",
    "      \"unit\": {\n",
    "        \"@id\": \"id.link/mm\"\n",
    "      }\n",
    "    },\n",
    "    {\n",
    "      \"@id\": \"d3382985\",\n",
    "      \"@type\": \"https://openminds.ebrains.eu/core/QuantitativeValue\",\n",
    "      \"value\": 0.485,\n",
    "      \"unit\": {\n",
    "        \"@id\": \"id.link/mm\"\n",
    "      }\n",
    "    },\n",
    "    {\n",
    "      \"@id\": \"84fc5d5d\",\n",
    "      \"@type\": \"https://openminds.ebrains.eu/core/QuantitativeValue\",\n",
    "      \"value\": 40.4396,\n",
    "      \"unit\": {\n",
    "        \"@id\": \"id.link/mm\"\n",
    "      }\n",
    "    }\n",
    "  ]\n",
    "}\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25abaf8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "view = plotting.plot_img(bigbrain_template.fetch(), cut_coords = tuple(point), cmap='gray')\n",
    "view.add_markers([tuple(point)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef21c2c3",
   "metadata": {},
   "source": [
    "## Probabilistic assignment of points\n",
    "\n",
    "Parcellation maps in siibra can assign their structures to locations and images.\n",
    "\n",
    "Note that siibra can warp locations between different template spaces. It does this automatically in many cases. For example, the Julich-Brain probability maps are defined in MNI space, while our point above is from BigBrain. siibra will detect this and convert accordingly for the assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "063cb42f",
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_pmaps.assign(point)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66e5057c",
   "metadata": {},
   "source": [
    "## Your turn!\n",
    "\n",
    "<div class='alert alert-block alert-info'>\n",
    "Pick a location from the online viewer, generate a siibra.Point and assign it to brain structures.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "024a800e",
   "metadata": {},
   "source": [
    "## Probabilistic assignment of image signals\n",
    "\n",
    "Just like points, siibra can assign brain regions to structures in an image volume. Here we load an example volume, and feed it to the same assign() method as the point. Since we have here image information, the asignment produces some scores for each match. siibra will split the volume into disconnected components. This typically results in many more assignments, so we will filter the resulting list using the scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e6908b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = siibra.volumes.from_file(\n",
    "    'ohbm-2023-example-input.nii.gz',\n",
    "    space='mni152',\n",
    "    name='Example input'\n",
    ")\n",
    "plotting.plot_stat_map(img.fetch())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f41f6b9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "assignments = julich_pmaps.assign(img)\n",
    "assignments[assignments.correlation > 0.3].sort_values(\"input containedness\")  # refer to pandas documentation for more..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17cc32d8",
   "metadata": {},
   "source": [
    "# Part IV - Collecting multimodal features\n",
    "\n",
    "\n",
    "`siibra` provides access to data features of different modalities. The features and their query functions are bundled in the module `siibra.features`. We can choose different types of features from this module. The feature types are organized in a hierarchy under the most abstract type `siibra.features.Feature`. All other feature types are subclasses of it. Let's look at this hierarchy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39f7cdc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "siibra.features.render_ascii_tree(\"Feature\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcc2d87e",
   "metadata": {},
   "source": [
    "### Densities of neurotransmitter receptors\n",
    "\n",
    "Features can be queried for brain regions, parcellations and location objects (such as volumes of interest in a refrence space) with `siibra.features.get()`, which accepts a query object and a feature type. It will query all subclasses of the given feature type, if any. Here is a simple example for getting a receptor density fingerprint in region V1. The data is of tabular type, provided as a pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e6834c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 = siibra.get_region(parcellation='julich 2.9', region='v1 left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc0bc915",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = siibra.features.get(\n",
    "    v1, siibra.features.molecular.ReceptorDensityFingerprint\n",
    ")\n",
    "\n",
    "# fetch the first one\n",
    "f = features[0]\n",
    "print(f.name)\n",
    "f.data.T  # we transpose the table for display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a07ded96",
   "metadata": {},
   "outputs": [],
   "source": [
    "# besides the data table, the features have some additional metadata\n",
    "print(f.modality)\n",
    "print(f.description)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0fa5b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62400435",
   "metadata": {},
   "source": [
    "### Understanding what was matched"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5998ce64",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.last_match_description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f918f06",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(julich_brain.get_region(\"occipital cortex\").tree2str())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74673b64",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = siibra.features.get(\n",
    "    julich_brain.get_region(\"occipital cortex\"),\n",
    "    siibra.features.molecular.ReceptorDensityFingerprint\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6fe4f47",
   "metadata": {},
   "outputs": [],
   "source": [
    "for f in features:\n",
    "    print(f.last_match_description)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf4fb5b2",
   "metadata": {},
   "source": [
    "## Your turn! \n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    Query layerwise cell densities for V1\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cfc44c9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "edfc1517",
   "metadata": {},
   "source": [
    "### Gene Expressions from the Allen Atlas\n",
    "\n",
    "siibra also implements a live queryto gene expression data from the Allen atlas to extract regional gene expression levels. Gene expressions are linked to atlas regions by coordinates of their probes in MNI space. When called with a brain region, siibra.features.get will generate a mask of this region in MNI space to filter the probes. It provides the regional expression levels as tabular data, together with the MNI coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d816d5c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = siibra.features.get(\n",
    "    v1, siibra.features.molecular.GeneExpressions, gene=[\"TAC1\", \"MAOA\", \"GABARAPL2\"]\n",
    ")\n",
    "fig = features[0].plot()\n",
    "features[0].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff6bb587",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can plot the MNI coordinates to localize the measures and verify they are located in V1.\n",
    "locations = siibra.PointSet(features[0].data.mni_xyz.tolist(), space=\"mni152\")\n",
    "\n",
    "from nilearn import plotting\n",
    "mask = v1.get_regional_map(\"mni152\")\n",
    "display = plotting.plot_glass_brain(mask.fetch(), cmap='viridis')\n",
    "display.add_markers(locations.as_list(), marker_size=2) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af798797",
   "metadata": {},
   "source": [
    "### Structural and functional connectivity\n",
    "\n",
    "`siibra` provides connectivity matrices with parcellation averaged structural and functional measurments for different subjects of different cohorts. Here we request some streamline counts for Julich Brain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a85c404",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = siibra.features.get(\n",
    "    siibra.parcellations.get('julich 2.9'),\n",
    "    siibra.features.connectivity.StreamlineCounts\n",
    ")\n",
    "f = features[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb2423ea",
   "metadata": {},
   "source": [
    "Let's check the cohort and subject ids of the first connectivity feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "349fbb43",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f.cohort)\n",
    "print(f.indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b150f07",
   "metadata": {},
   "source": [
    "We can retrieve the matrix of a single subject using `get_matrix()`. If we leave the subject specification out, `siibra` will compute the mean matrix across subjects.\n",
    "\n",
    "Again, the result is a pandas dataframe, with the notable property that the row and column indices are full region objects for further reference. This implies in particular, that we can directly associate each measure with the corresponding information in the parcellation, and with a mask of a parcellation map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dce2d15",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc_hcp_000 = f.get_element('000')\n",
    "sc_hcp_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9da76b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also plot the matrix\n",
    "sc_hcp_000.plot(logscale=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b1f53b0",
   "metadata": {},
   "source": [
    "As an example, we retrieve the centroids in MNI152 space and plot the connnectivity graph in 3D:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f421e7fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "node_coords = sc_hcp_000.compute_centroids(space='mni152')\n",
    "plotting.view_connectome(\n",
    "    adjacency_matrix=sc_hcp_000.data,\n",
    "    node_coords=node_coords,\n",
    "    edge_threshold=\"99%\",\n",
    "    node_size=3, colorbar=False,\n",
    "    edge_cmap=\"bwr\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bacdf9c",
   "metadata": {},
   "source": [
    "## Working with BigBrain data\n",
    "\n",
    "### Extracting chunks from the 20 micron model\n",
    "\n",
    "To access full resolution data of BigBrain, we specify a bounding box in the physical space. For now, we just define a volume of interest from two corner points in the histological space. We specify the points with a string representation, which could be conveniently copy pasted from the interactive viewer siibra explorer. Note that the coordinates can be specified by 3-tuples, and in other ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4814d107",
   "metadata": {},
   "outputs": [],
   "source": [
    "voi = siibra.locations.BoundingBox(\n",
    "    point1=\"-30.590mm, 3.270mm, 47.814mm\",\n",
    "    point2=\"-26.557mm, 6.277mm, 50.631mm\",\n",
    "    space='bigbrain'\n",
    ")\n",
    "\n",
    "# Note: we can reuse the template object from above,\n",
    "# just fetch with different parameters\n",
    "bigbrain_chunk = bigbrain_template.fetch(voi=voi, resolution_mm=0.02)\n",
    "plotting.plot_img(bigbrain_chunk, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b85aa570",
   "metadata": {},
   "source": [
    "### Loading cortical profiles (by Wagstyl et al.)\n",
    "\n",
    "Cortical staining profiles in BigBrain have been precomputed by Konrad Wagstyl and colleagues. Siibra has integrated those as a feature type, for the left hemishere. We can thus query profiles by region!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "938b61ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = siibra.features.get(\n",
    "    siibra.get_region('julich', '4p left'),\n",
    "    siibra.features.cellular.BigBrainIntensityProfile\n",
    ")\n",
    "print(len(features))\n",
    "f = features[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81a6c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f140261d",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1eacea",
   "metadata": {},
   "source": [
    "### 1 micron sections\n",
    "\n",
    "HIBALL has released a range of 1 micron scans of BigBrain sections across the brain. Siibra can find those as VolumeOfInterest features. The result is a high-resolution image structure, just like the bigbrain template."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5020f4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "hoc5l = siibra.get_region('julich', 'hoc5 left')\n",
    "features = siibra.features.get(\n",
    "    hoc5l,\n",
    "    siibra.features.cellular.CellbodyStainedSection\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56655de6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's see the names of the found features\n",
    "for f in features:\n",
    "    print(f.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab643b0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's fetch the 1 micron section at a lower resolutioin, and display in 3D space.\n",
    "section1402 = features[3]\n",
    "plotting.plot_img(\n",
    "    section1402.fetch(),\n",
    "    bg_img=bigbrain_template.fetch(),\n",
    "    title=\"#1402\",\n",
    "    cmap='gray'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f97dbec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's fetch a crop inside hoc5 at full resolution.\n",
    "# we intersect the bounding box of hoc5l and the section\n",
    "hoc5_bbox = hoc5l.get_boundingbox('bigbrain').intersection(section1402)\n",
    "print(f\"Size of the bounding box: {hoc5_bbox.shape}\")\n",
    "\n",
    "# this is quite large, so we shrink it\n",
    "voi = hoc5_bbox.zoom(0.1)\n",
    "crop = section1402.fetch(voi=voi, resolution_mm=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab22d3e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_img(crop, bg_img=None, cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e38eb5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_img(crop, bg_img=bigbrain_template.fetch(), cmap='magma')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f5ac72f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}