neurolang_tutorial_01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # NeuroLang is a logic language enabling rich meta-analyses.\n",
    "\n",
    "The focus of neurolang is to express complex hypotheses and fit them to the data in a tool that is closer to human discourse. This tutorial is focused on the logic language aspect (more natural inputs coming!!). Neurolang is mainly aimed at the probabilistic and noisy aspect of relating neuroimaging information with well-formulated hypotheses that are expressible in natural languages.\n",
    "\n",
    "\n",
    " For instance, in the following image, a meta-analysis of the Neurosynth database was able to pinpoint the\n",
    " most likely cognitive terms associated with different functional networks (Abdallah et al [Sci Rep 2022](https://www.nature.com/articles/s41598-022-21801-4), [eLife 2022](https://elifesciences.org/articles/76926)).\n",
    "\n",
    " <img src=\"abdallah_et_al_fig2.png\" width=\"750\"/>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "assert(sys.version_info.major == 3 and 7 < sys.version_info.minor < 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install -q \"matplotlib<3.7\"\n",
    "%pip install -q git+https://github.com/wannesm/PySDD.git#egg=PySDD\n",
    "%pip install -q neurolang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "from neurolang.frontend import NeurolangPDL\n",
    "from neurolang.utils.server.engines import NeurosynthEngineConf\n",
    "import os\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "from typing import Callable\n",
    "from nilearn import plotting, datasets\n",
    "from matplotlib import pyplot as plt\n",
    "import nibabel as nib\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # Let's learn neurolang by extracting information from NeuroSynth!\n",
    "\n",
    "\n",
    " Let's start by creating NeuroLang instance. Specifically one with the Neurosynth database loaded. This database has three tables:\n",
    "\n",
    " * `TermInStudyTFIDF(term, tfidf, study)` which contains all terms mentioned in a study with their tfidf\n",
    " * `PeaksReported(x, y, z, study)` which contains the peaks, in MNI 152 space, reported in a study.\n",
    " * `Voxels(x, y, z)` which contains all voxels of the MNI 152 space, at 1mm resolution\n",
    " * `Study(s)` the set of all studies\n",
    " * `SelectedStudy(s)` which represents that every study in the database is an equiprobable token of neuroscientific knowledge. Specifically, the table represents the probabilistic events \"one study `s` was selected\" where each event has the probability `1 / <total number of Studies>` and these events are exclusive.\n",
    "\n",
    " The set of tables in a NeuroLang instance can be called the _extensional database_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neurolang_engine_factory = NeurosynthEngineConf(Path(os.getcwd()))\n",
    "nl = neurolang_engine_factory.create()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Let's make a first query obtaining all terms mentioned in all studies. Specifically we are going to encode the rule:\n",
    "\n",
    " Every `term` in the set `Term` is the first element of the tuples in the set `TermInStudyTFIDF`\n",
    "\n",
    "\n",
    " $\\forall term: Term(term) \\leftarrow TermInStudyTFIDF(term, tfidf, study)$\n",
    "\n",
    " and then we will add a line expressing that the expected answer is all the terms in the set `Term`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        Term(term) :- TermInStudyTFIDF(term, tfidf, study)\n",
    "        ans(term) :- Term(term)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Next, we will extract the top 99.999% most important terms according to their TFIDF. For this we will add a `numpy` function into our engine.\n",
    " Then, we will define the TFIDF threshold, and use it to select the terms with two rules representing\n",
    "\n",
    " \"Define all Important Terms as those whose TFIDF is in the top 99.999%\":\n",
    "\n",
    " * The TFIDFThreshold is the TFIDF valuee at the 99.999 percentile\n",
    " * The ImportanTerms are those such that there is at least one article where the TFIDF is larger than the threshold\n",
    "   * $\\forall term, tfidf: ImportantTerm(term, tfidf)\\leftarrow \\exists study: TermInStudy(term, tfidf, study) \\wedge tfidf > thr \\wedge TFIDFThreshold(thr)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nl.add_symbol(\n",
    "    np.percentile, name=\"percentile\", type_=Callable[[float, float], float]\n",
    ")\n",
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        TFIDFThreshold(percentile(tfidf, 99.999)) :- TermInStudyTFIDF(t, tfidf, study)\n",
    "        ImportantTerm(term, tfidf) :- TermInStudyTFIDF(term, tfidf, study) & (tfidf > thr) & TFIDFThreshold(thr)\n",
    "        ans(term, tfidf) :- ImportantTerm(term, tfidf)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res.as_pandas_dataframe().sort_values(by='tfidf'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " We can do the same thing but using more convenient expressions in the query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        Threshold_tfidf(percentile(tfidf, (100 - .001))) :- TermInStudyTFIDF(t, tfidf, study)\n",
    "        Term(t, tfidf) :- TermInStudyTFIDF(t, tfidf, study) & (tfidf > thr) & Threshold_tfidf(thr)\n",
    "        ans(term, tfidf) :- Term(term, tfidf)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res.as_pandas_dataframe().sort_values(by='tfidf'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Next we are going to extract all MNI voxels that have been reported active in a study mentioning the term \"emotion\".\n",
    "\n",
    " Specifically,\n",
    " 1. we will consider active those who are at most within 4mm of a peak reported in a study.\n",
    "      * \"A Voxel is Reported in a study if there is a peak reported in such study, and the voxel is within 4mm of this peak\n",
    "\n",
    "\n",
    " 2. we will filter all the studies according to the importance of the word \"emotion\" in the study\n",
    "      * \"a term is in a emotion study if the therm is mentioned in the study and it's TFIDF in the study is larger tham 0.1\"\n",
    "\n",
    "\n",
    " 3. We will count how many times a voxel was reported in a emotion study\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        VoxelReported (x, y, z, study) :- PeakReported(x2, y2, z2, study) & Voxel(x, y, z) & (d == EUCLIDEAN(x, y, z, x2, y2, z2)) & (d < 4)\n",
    "        TermInEmotionStudy(term, study) :- TermInStudyTFIDF(term, tfidf, study) & (tfidf > 0.1)\n",
    "        ReportedInEmotionStudy(x, y, z, count(study)) :- VoxelReported (x, y, z, study) & TermInEmotionStudy(\"emotion\", study)\n",
    "        ans(x, y, z, s) :- ReportedInEmotionStudy(x, y, z, s)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # Creating overlays as results for the analysis\n",
    "\n",
    " The results as a voxel list are not very easy to understand. So we will agregate all voxel counts in an image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        VoxelReported (x, y, z, study) :- PeakReported(x2, y2, z2, study) & Voxel(x, y, z) & (d == EUCLIDEAN(x, y, z, x2, y2, z2)) & (d < 4)\n",
    "        TermInEmotionStudy(term, study) :- TermInStudyTFIDF(term, tfidf, study) & (tfidf > 0.1)\n",
    "        ReportedInEmotionStudy(x, y, z, count(study)) :- VoxelReported (x, y, z, study) & TermInEmotionStudy(\"emotion\", study)\n",
    "        EmotionCountImage(agg_create_region_overlay(x, y, z, s)) :- ReportedInEmotionStudy(x, y, z, s)\n",
    "        ans(img) :- EmotionCountImage(img)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res)\n",
    "region = res.as_pandas_dataframe().iloc[0, 0]\n",
    "thr = np.percentile(region.overlay, 95)\n",
    "plotting.plot_img(\n",
    "    region.spatial_image(), colorbar=True,\n",
    "    threshold=thr, bg_img=datasets.load_mni152_template(),\n",
    "    display_mode='y', cut_coords=np.linspace(-50, 50, 10)\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Excercise: repeat the task but counting the a study reports a voxel when the study mentions \"language\" and \"audition\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # Let's now talk about probabilities! What is the probability that a voxel is mentioned in a Study?\n",
    "\n",
    " Now let's formulate the same question but in a probabilistic manner. Specifically we will produce the probability that a voxel is mentioned in a study, given that the study mentions the word \"emotion\", in shorthand notation this is\n",
    " * $P(voxel \\in Study | ``emotion\" \\in Study)$\n",
    " Which we can write in neurolang as:\n",
    " * `VoxelsGivenStudyMentionsEmotion(x, y, z, PROB) :- VoxelReported(x, y, z, study) // (TermInEmotionStudy(\"emotion\", study) & SelectedStudy(study))`\n",
    "   In this case the `PROB` term in the left-hand side of the rule, or head, means that we will obtain the probability of the stochastic event on the right. The `//` operator means conditional probability, and `SelectedStudy` as we described at the beginning of the notebook, is a probabilistic event, denoting that the study `study` is mentioned with probability `1 / <total number of studies>.`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        VoxelReported (x, y, z, study) :- PeakReported(x2, y2, z2, study) & Voxel(x, y, z) & (d == EUCLIDEAN(x, y, z, x2, y2, z2)) & (d < 4)\n",
    "        TermInEmotionStudy(term, study) :- TermInStudyTFIDF(term, tfidf, study) & (tfidf > 0.1)\n",
    "        VoxelsGivenStudyMentionsEmotion(x, y, z, PROB) :- VoxelReported(x, y, z, study) // (TermInEmotionStudy(\"emotion\", study) & SelectedStudy(study))\n",
    "        image(agg_create_region_overlay(x, y, z, p)) :- VoxelsGivenStudyMentionsEmotion(x, y, z, p)\n",
    "        ans(img) :- image(img)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res)\n",
    "region = res.as_pandas_dataframe().iloc[0, 0]\n",
    "thr = np.percentile(region.overlay, 95)\n",
    "plotting.plot_stat_map(region.spatial_image(), threshold=thr, display_mode='y', cut_coords=np.linspace(-50, 50, 10))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some hints:\n",
    "\n",
    "In [V. Iovene's thesis](https://theses.hal.science/tel-03653150) you can find the implementation of different Coordinate-Based Meta-Analyses analyses in NeuroLang."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Excercise:\n",
    "\n",
    " * Produce the image of the probability that a voxel will be mentioned \"language\" and \"audition\" being mentioned simultaneously on a study\n",
    " * Produce the image of the log likelihood ratios comparing the probability that a voxel will be mentioned in a study mentioning \"language\" and that it will be reported in a study mentioning \"audition\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # Reverse Inference. Figuring out the cognitive terms associated to the Yeo 7 Network parcellation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make the atlas into a set of tuples (x, y, z, label)\n",
    "yeo_atlas = datasets.fetch_atlas_yeo_2011()\n",
    "yeo7 = nib.load(yeo_atlas['thin_7'])\n",
    "plotting.plot_roi(yeo_atlas['thin_7'], colorbar=True)\n",
    "\n",
    "yeo_data = yeo7.get_fdata()\n",
    "yeo_ijk = np.nonzero(yeo_data)\n",
    "yeo_xyz = nib.affines.apply_affine(yeo7.affine, np.transpose(yeo_ijk)[:, :3])\n",
    "yeo_labels = yeo_data[yeo_ijk]\n",
    "yeo7_melt = np.c_[yeo_xyz, yeo_labels]\n",
    "\n",
    "# Add the Set to the neurolang engine as the \"Yeo7\" set\n",
    "\n",
    "nl.add_tuple_set(yeo7_melt, name='Yeo7')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        NetworkReported(network, study) :- Yeo7(x2, y2, z2, network) & PeakReported(x, y, z, study) & (d == EUCLIDEAN(x, y, z, x2, y2, z2)) & (d <= 2)\n",
    "        Network(network) :- Yeo7(..., ..., ..., network) \n",
    "        TermInStudy(term, study) :- TermInStudyTFIDF(term, tfidf, study) & (tfidf > 0.1)\n",
    "        NetworkTerm(term, network, PROB) :- TermInStudy(term, study) // (NetworkReported(network, study) & SelectedStudy(study))\n",
    "        Term(term, PROB) :- TermInStudy(term, study) & SelectedStudy(study)\n",
    "        NetworkLR(term, network, ll) :- NetworkTerm(term, network, p) & Term(trem, p1) & (ll == log(p / p1))\n",
    "        ans(term, network, p) :- NetworkLR(term, network, p)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res.as_pandas_dataframe().sort_values(by='p').groupby('network').tail(2).sort_values(\"network\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with nl.scope as e:\n",
    "    res = nl.execute_datalog_program(\n",
    "    \"\"\"\n",
    "        NetworkReported(network, study) :- Yeo7(x2, y2, z2, network) & PeakReported(x, y, z, study) & (d == EUCLIDEAN(x, y, z, x2, y2, z2)) & (d <= 2)\n",
    "        Network(network) :- Yeo7(..., ..., ..., network) \n",
    "        NetworkExclusiveReported(network, study) :- NetworkReported(network, study) & ~exists(n2; Network(n2) & (n2 != network) & NetworkReported(n2, study))\n",
    "        TermInStudy(term, study) :- TermInStudyTFIDF(term, tfidf, study) & (tfidf > 0.1)\n",
    "        NetworkTerm(term, network, PROB) :- TermInStudy(term, study) // (NetworkExclusiveReported(network, study) & SelectedStudy(study))\n",
    "        Term(term, PROB) :- TermInStudy(term, study) & SelectedStudy(study)\n",
    "        NetworkLR(term, network, ll) :- NetworkTerm(term, network, p) & Term(trem, p1) & (ll == log(p / p1))\n",
    "        ans(term, network, p) :- NetworkLR(term, network, p)\n",
    "    \"\"\"\n",
    "    )\n",
    "    display(res.as_pandas_dataframe().sort_values(by='p').groupby('network').tail(2).sort_values(\"network\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Excercise: Do the same but with the Yeo 17 network set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yeo_atlas = datasets.fetch_atlas_yeo_2011()\n",
    "yeo17 = nib.load(yeo_atlas['thin_17'])\n",
    "plotting.plot_roi(yeo_atlas['thin_17'], colorbar=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "neurolang_tutorial_310a",
   "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.10.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}