notebooks/13-diskussion-und-ausblick.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ein paar interessante Studien"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### module importieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn import preprocessing\n",
    "\n",
    "from nilearn import plotting, image, input_data\n",
    "from nistats import first_level_model, reporting, regression\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set_context('talk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Haxby et al. 2001 Science\n",
    "\n",
    "Dekodieren von Hirnaktivität mittels Ähnlichkeitsmaßen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(\"https://science.sciencemag.org/content/sci/293/5539/2425.full.pdf\", width=900, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### unsere daten laden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_img_list = glob.glob('../data/sub-01_task-imagery_run-0*_bold_space-MNI152NLin2009cAsym_preproc.nii.gz')\n",
    "fmri_img_list.sort()\n",
    "fmri_img_list = fmri_img_list[1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### extrahieren der Daten mit einem Masker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_img = image.mean_img([image.mean_img(f) for f in fmri_img_list])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masker = input_data.NiftiMasker(standardize=True,\n",
    "                                detrend=True,\n",
    "                                smoothing_fwhm=6,\n",
    "                                mask_strategy='epi')\n",
    "masker.fit(mean_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_roi(masker.mask_img_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### jeden Run einzeln extrahieren, dann an die anderen anhängen\n",
    "\n",
    "So erreichen wir, dass jeder Run einzeln z-transformiert wird. Dadurch ist auch der Mittlwert jedes Runs bei null und wir brauchen keine Konstante mehr."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_big_df(fmri_img_list,masker):\n",
    "    big_df = pd.DataFrame()\n",
    "    for im in fmri_img_list:\n",
    "        data = masker.transform(im)\n",
    "        df = pd.DataFrame(data)\n",
    "        big_df = pd.concat([big_df,df])\n",
    "    return big_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_df = make_big_df(fmri_img_list,masker)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### design laden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "onset_csvs = glob.glob('../data/sub-01_task-imagery_run-0*_events.tsv')\n",
    "onset_csvs.sort()\n",
    "onset_csvs = onset_csvs[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_onsets(csv_file):\n",
    "    \n",
    "    df = pd.read_csv(csv_file,sep='\\t')\n",
    "    df.loc[:,'trial_num'] = [('0%s'%x)[-2:] for x in df.index.get_values().astype(str)]\n",
    "    df.loc[:,'trial_type'] = df.loc[:,'trial_domain']\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_volumes = image.load_img(fmri_img_list[-1]).shape[-1]\n",
    "my_tr = 3\n",
    "my_duration = n_volumes*my_tr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_design(df,my_tr=my_tr,my_duration=my_duration):\n",
    "    frame_times = np.arange(0,my_duration,my_tr)\n",
    "    design_matrix = first_level_model.make_first_level_design_matrix(frame_times=frame_times,\n",
    "                                                                     events=df,\n",
    "                                                                     drift_model=None\n",
    "                                                                    )\n",
    "    return design_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_big_design(onset_csvs):\n",
    "    \n",
    "    design_matrix = pd.DataFrame()\n",
    "    for run_num,onset_csv in enumerate(onset_csvs):\n",
    "        this_onset = make_onsets(onset_csv)\n",
    "        this_design = make_design(this_onset).iloc[:,:-1] # konstante weg\n",
    "        this_design.index = this_design.index+my_duration*run_num\n",
    "        this_design.columns = ['%s-run%s'%(c,run_num+1) for c in this_design.columns]\n",
    "\n",
    "        design_matrix = pd.concat([design_matrix, this_design],sort=False)\n",
    "        design_matrix = design_matrix.fillna(0)\n",
    "\n",
    "    return design_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "design_matrix = make_big_design(onset_csvs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1,figsize=(14,10))\n",
    "sns.heatmap(design_matrix,ax=ax,cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Daten und Design\n",
    "\n",
    "Hier multiplizieren wir den Zeitverlauf jedes Voxels mit der erwarteten HRF jeder Bedingung; so bekommen wir eine Art gewichteten Mittelwert für jede Bedingung."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "this_design = design_matrix.iloc[:200,4].values\n",
    "this_voxel = big_df.iloc[:200,4000].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,6))\n",
    "plt.plot(this_design,label='model of hrf for one condition')\n",
    "plt.plot(this_voxel,label='timecourse in one voxel')\n",
    "plt.xlabel('time')\n",
    "plt.legend()\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,6))\n",
    "plt.plot(this_design * this_voxel,label='timecourse in one voxel for one condition')\n",
    "plt.xlabel('time')\n",
    "plt.legend(loc=(0.5,1))\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### gemittelte Aktivität für jede Bedingung"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def make_block_df(big_df,design_matrix):\n",
    "    block_df = pd.DataFrame()\n",
    "    for cond in design_matrix.columns:\n",
    "        print(cond)\n",
    "        f = lambda x:x*design_matrix[cond].values\n",
    "        mean_block = big_df.apply(f).mean()\n",
    "        mean_block_df = pd.DataFrame(mean_block).T\n",
    "        mean_block_df.index = [cond]\n",
    "        block_df = pd.concat([block_df, mean_block_df])\n",
    "    return block_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "block_df = make_block_df(big_df,design_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "block_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "block_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beispiel: Die Daten aus einer Zeile in den Hirnraum zurückführen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "block_im = masker.inverse_transform(block_df.loc['SPRACHE-run3',:])\n",
    "plotting.view_img(block_im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aufteilen in zwei Hälften"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_half = block_df.iloc[:10,:]\n",
    "first_half"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_half = block_df.iloc[10:,:]\n",
    "second_half"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### erste Hälfte als Modell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_half.index = pd.MultiIndex.from_tuples([x.split('-') for x in first_half.index])\n",
    "first_half"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_half = first_half.groupby(level=0).mean()\n",
    "first_half"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alle Blocks der 2. Hälfte mit denen der 1. Hälfte Korrelieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_corr_df = pd.DataFrame()\n",
    "for i in first_half.index:\n",
    "    this_corr = second_half.T.corrwith(first_half.T.loc[:,i])\n",
    "    this_corr_df = pd.DataFrame(this_corr).T\n",
    "    this_corr_df.index = [i]\n",
    "    big_corr_df = pd.concat([big_corr_df,this_corr_df])\n",
    "big_corr_df = big_corr_df.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1,figsize=(10,6))\n",
    "sns.heatmap(big_corr_df,annot=True,ax=ax,cmap='RdBu_r',vmin=-1,vmax=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kriegeskorte et al. 2008 Neuron\n",
    "\n",
    "Representational Similiarity Analysis (RSA). Wie ähnlich ist die \"Repräsentation\" verschiedener kognitiver Zustände im Gehirn bzw. bestimmten Hirnregionen?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IFrame(\"https://www.cns.nyu.edu/kianilab/papers/Kriegeskorte_Kiani_Bandettini.pdf\", width=900, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dissimilarity berechnen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_df = block_df.T.corr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(corr_df,cmap='RdBu_r',vmin=-1,vmax=1,square=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diss_df = 1 - corr_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(diss_df,cmap='Reds',vmin=0,vmax=2,square=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.clustermap(diss_df,figsize=(16,16),cmap='Reds',annot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-Dimensionale Skalierung"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import manifold\n",
    "from sklearn.metrics import euclidean_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds = manifold.MDS(n_components=2,\n",
    "                   dissimilarity=\"precomputed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds_positions = mds.fit( diss_df ).embedding_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds_df = pd.DataFrame(mds_positions,index=diss_df.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mds_df.loc[:,'group'] = [x.split('-')[0] for x in mds_df.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1,figsize=(12,12))\n",
    "sns.scatterplot(data=mds_df,x=0,y=1,hue='group',ax=ax)\n",
    "for i in mds_df.index:\n",
    "    x,y = mds_df.loc[i,[0,1]]\n",
    "    plt.annotate(s=i,xy=(x,y))\n",
    "#plt.xticks([]); plt.yticks([])\n",
    "#plt.xlabel(''); plt.ylabel('')\n",
    "#sns.despine(left=True,bottom=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hasson et al. 2004 Science\n",
    "\n",
    "Reverse Correlation: Welche Art von Stimulus aktiviert eine Hirnregion maximal?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IFrame(\"http://www.psiexp.ss.uci.edu/research/teachingP140C/Papers/hasson_2004.pdf\", width=900, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Daten in Hirnraum zurückführen und Signal aus Koordinate extrahieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_im = masker.inverse_transform(big_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere_masker = input_data.NiftiSpheresMasker([(-52,12,28)],\n",
    "                                              radius=5,\n",
    "                                              t_r=my_tr,\n",
    "                                              low_pass=0.1).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere_data = sphere_masker.transform(big_im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere_df = pd.DataFrame(sphere_data)\n",
    "sphere_df.columns = [','.join([str(i) for i in x]) for x in sphere_masker.seeds_]\n",
    "sphere_df.loc[:,'time'] = sphere_df.index*my_tr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere_df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_name = sphere_df.columns[0]\n",
    "col_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aktivität in der extrahieren Region über die Zeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1,figsize=(16,6))\n",
    "ax.plot(sphere_df.loc[:,col_name])\n",
    "plt.axhline(0,color='k')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('z-transformed signal')\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Höchste Werte"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "my_max_values = sphere_df.loc[:,col_name].sort_values(ascending=False)[:10].values\n",
    "my_max_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Index der höchsten Werte"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_max_idx = sphere_df.loc[:,col_name].sort_values(ascending=False)[:10].index\n",
    "my_max_idx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zeit in Sekunden der höchsten Werte"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_max_time = [ sphere_df.loc[i,'time'] for i in my_max_idx]\n",
    "my_max_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inhalt der Design-Matrix zu diesen Zeitpunkten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_max_label = [design_matrix.loc[i].idxmax() for i in my_max_time]\n",
    "my_max_label"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Annotierter Zeitverlauf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,1,figsize=(16,6))\n",
    "ax.plot(sphere_df.loc[:,col_name],alpha=0.9)\n",
    "for (s,x,y) in zip(my_max_label,my_max_idx,my_max_values):\n",
    "    plt.plot(x,y,'o',color='r',alpha=0.5)\n",
    "    plt.annotate(s=s,xy=(x,y))\n",
    "plt.axhline(0,color='k')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('z-transformed signal')\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DIY fMRI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IFrame('https://neuwritesd.org/2015/07/10/got-a-computer-we-have-data-lets-do-neuroscience/', width=700, height=350)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "IFrame('https://antmelder.files.wordpress.com/2014/06/tumblr_m4j2quesvw1qedj2ho1_1280.jpg', width=700, height=500)"
   ]
  }
 ],
 "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
}