notebooks/08-whole-brain-analyse-mit-hrf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Whole-Brain Ergebniskarten mit HRF als Modell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### module importieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from nilearn import image, input_data, plotting\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from scipy import stats\n",
    "from nistats import hemodynamic_models\n",
    "\n",
    "sns.set_context('poster')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dateiname des funktionellen Bildes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epi_file = '../data/sub-01_task-imagery_run-02_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extraktion der Daten aus der Region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whole-Brain Maske definieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_standard = True\n",
    "my_detrend = True\n",
    "my_smoo = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain_masker = input_data.NiftiMasker(mask_img='../data/sub-01_task-imagery_run-02_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz',\n",
    "                                            standardize=my_standard,\n",
    "                                            detrend=my_detrend,\n",
    "                                            smoothing_fwhm=my_smoo).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Daten auslesen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain_data = brain_masker.transform(epi_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Daten in Tabelle packen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain_df = pd.DataFrame(brain_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Studiendesign laden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "design_df = pd.read_csv('../data/sub-01_task-imagery_run-02_events.tsv',sep='\\t',index_col=[2,4])\n",
    "design_df = design_df.sort_index()\n",
    "design_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modell der hämodynamischen Antwort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_tr = 3.0\n",
    "my_onset = 0.0\n",
    "my_length = 32.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrf_model = hemodynamic_models.spm_hrf(tr=my_tr,\n",
    "                                   oversampling=my_tr,\n",
    "                                   time_length=my_length+my_onset,\n",
    "                                   onset=my_onset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(hrf_model,'-o',label='erwartete messbare Antwort')\n",
    "plt.axvline(my_onset,color='r',label='neuronales Feuern')\n",
    "plt.axhline(0,color='k')\n",
    "sns.despine()\n",
    "plt.legend(loc=(1,0.5))\n",
    "plt.xlabel('Zeit in Sekunden')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Design und HRF miteinander verbinden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_volumes = image.load_img(epi_file).shape[-1]\n",
    "number_of_volumes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_onsets(df,number_of_volumes=number_of_volumes,my_tr=my_tr):\n",
    "    design = np.zeros(int(number_of_volumes*my_tr))\n",
    "    # für jede Zeile in der Tabelle\n",
    "    for start in range(df.shape[0]):\n",
    "        # Inhalte der Zeile\n",
    "        block = df.iloc[start]\n",
    "        # Onsets und Durations extrahieren\n",
    "        block_onset = block['onset']\n",
    "        block_duration = block['duration']\n",
    "        # Ab Onset so viele 1en einfügen wie der Block lang ist\n",
    "        for event in range(block_duration):\n",
    "            design[block_onset+event] = 1.\n",
    "    return design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_design(my_condition):\n",
    "    \n",
    "    beispiel_design = make_onsets(design_df.loc[my_condition,['onset','duration']])\n",
    "    hrf_design_in_tr = np.convolve(beispiel_design,hrf_model)[::int(my_tr)]\n",
    "    \n",
    "    return hrf_design_in_tr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hrf_design_in_tr = my_design('MOTORIK')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hrf_design_in_tr);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Design für alle Bedingungen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Leere Tabelle in die wir hineinschreiben\n",
    "design_matrix_df = pd.DataFrame()\n",
    "# für jede der 5 Bedingungen aus design_df\n",
    "for condition in design_df.index.levels[0]:\n",
    "    # verbinde das design mit der hrf\n",
    "    hrf_design_in_tr = my_design(condition)\n",
    "    # mach daraus eine tabelle\n",
    "    hrf_df = pd.DataFrame(hrf_design_in_tr)\n",
    "    # der spaltenname der tabelle entspricht dem namen der bedingung\n",
    "    hrf_df.columns = [condition]\n",
    "    # die akuelle tabelle wird der großen tabelle angehängt\n",
    "    design_matrix_df = pd.concat([design_matrix_df,hrf_df],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "design_matrix_df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualisierung mit der Zeit auf der x-Achse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,6))\n",
    "for condition in design_matrix_df.columns:\n",
    "    plt.plot(design_matrix_df.loc[:,condition],label=condition)\n",
    "sns.despine(trim=True)\n",
    "plt.legend(loc=(0.95,0.3))\n",
    "plt.xlabel('time in volumes')\n",
    "plt.ylabel('predicted activity\\n(arbitrary units)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative Ansicht \"von oben\" als Heatmap (Zeit auf der y-Achse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1,figsize=(7,9))\n",
    "sns.heatmap(design_matrix_df,\n",
    "            cbar_kws={'label': 'predicted activity\\n(arbitrary units)'},\n",
    "            cmap='Greys_r',\n",
    "            ax=ax)\n",
    "plt.ylabel('time in volumes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Korrelationen rechnen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bedingung auswählen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_onsets = design_matrix_df.loc[:,'MOTORIK']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_onsets.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Korrelationen des Signalvarlaufs jedes Voxels mit der Bedingung rechnen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_data = brain_df.corrwith(my_onsets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nachdem der Zeitverlauf (203 Volumen) jedes Voxels mit dem HRF-Modell korreliert wurde, gibt es so viele Korrelationen wie es Voxel gibt:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "corr_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Korrelationen in Hirnbild umwandeln (den Korrelationswert jedes Voxels an die entsprechende Stelle im 3D-Hirnraum packen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_img = brain_masker.inverse_transform(corr_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(corr_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hirnbild visualisieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(corr_img,\n",
    "                       threshold=0.3,\n",
    "                       display_mode='z',\n",
    "                       cut_coords=8,\n",
    "                       bg_img='../data/sub-01_T1w_space-MNI152NLin2009cAsym_preproc.nii.gz');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (seminar)",
   "language": "python",
   "name": "seminar"
  },
  "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}