5_tSNR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In this notebook, we'll convert the raw timeseries to tSNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import glob\n",
    "import nibabel as nb\n",
    "from scipy.stats import stats\n",
    "from tqdm import tqdm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define where wb_command is installed:\n",
    "wb_home = '/home/vanes/Downloads/workbench/bin_rh_linux64/'\n",
    "ts_home_dir = '/home/shared/2018/visual/HCP7TFIXED/'\n",
    "\n",
    "data_dir = '/home/shared/2018/visual/hcp_cerebellum/'\n",
    "ts_write_dir = os.path.join(data_dir,'timeseries')\n",
    "resource_dir = os.path.join('/home','vanes','git','hcp_cerebellum','resources')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ts_files = glob.glob(os.path.join(ts_home_dir,'*'))\n",
    "ss=os.listdir(ts_home_dir)\n",
    "sd=np.array([os.path.isdir(os.path.join(ts_home_dir,s)) for s in ss])\n",
    "dirs = np.array(ss)[sd]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loop over subjects and save subcortical niftis\n",
    "# generic wb command:\n",
    "wbc = wb_home + 'wb_command -cifti-separate {cii} COLUMN -volume-all {cii_n}.nii.gz'\n",
    "\n",
    "for si, s in enumerate(np.sort(dirs)):\n",
    "\n",
    "    fns = glob.glob(os.path.join(ts_home_dir,s,'*_clean.dtseries.nii'))\n",
    "    for fn in fns:\n",
    "        task = fn.split('/')[-1].split('_')[1]\n",
    "        # convert to subcortical nifti\n",
    "        save_fn =os.path.join(ts_write_dir,'sub_%d_%s'%(si,task))\n",
    "        wbc_c = wbc.format(cii=fn, cii_n=save_fn)\n",
    "        os.system(wbc_c)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dims = {\n",
    "    'ang':0,\n",
    "    'ecc':1,\n",
    "    'gain':2,\n",
    "    'meanvol':3,\n",
    "    'r2':4,\n",
    "    'rfsize':5\n",
    "}\n",
    "\n",
    "def create_roi_mask():\n",
    "    \"\"\"\n",
    "    This mask deselects all voxels outside the retinotopic\n",
    "    clusters identified in the average subject.\n",
    "    \"\"\"\n",
    "    \n",
    "    roi_fn = os.path.join(resource_dir,'volume_masks','cerebellum_retmaps.nii')\n",
    "    maskimg = nb.load(roi_fn)\n",
    "    roimask = (maskimg.get_data()>0) # valid voxels  \n",
    "\n",
    "    return roimask\n",
    "\n",
    "def determine_best_subjects():\n",
    "    \n",
    "    \"\"\"\n",
    "    This function determines the best subjects based on the median\n",
    "    r-squared within the retinotopic clusters defined in the average subject.\n",
    "    \"\"\"\n",
    "\n",
    "    roimask = create_roi_mask()\n",
    "\n",
    "    all_r2 = []\n",
    "    for sj in range(181):      \n",
    "\n",
    "        # load the prf results nifti\n",
    "        fn = os.path.join(data_dir,'all_subjects','prfresults_subject_%d.dscalar_data_sub.nii.gz'%sj)\n",
    "        img = nb.load(fn)\n",
    "        data = img.get_data()\n",
    "\n",
    "        # mask r2s with roi mask:\n",
    "        r2s = np.ravel(data[roimask,dims['r2']])    \n",
    "        # and get median within these voxels\n",
    "        all_r2.append(np.nanmedian(r2s))\n",
    "\n",
    "    # sort subjects based on median r2\n",
    "    best_subjects = np.argsort(all_r2)[::-1]\n",
    "\n",
    "    return best_subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s = determine_best_subjects()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/vanes/anaconda2/envs/analysis/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in divide\n"
     ]
    }
   ],
   "source": [
    "tasks = ['RETCON','RETEXP','RETBAR1','RETBAR2','RETCW','RETCCW']\n",
    "# loop over subjects and compute tSNR\n",
    "for si, sj in zip(np.arange(181),np.array(s)[these_s]):\n",
    "    tsnrs = []\n",
    "    for task in tasks:\n",
    "        fn = os.path.join(ts_write_dir,'sub_%d'%sj+'_'+task+'.nii.gz')\n",
    "        img = nb.load(fn)\n",
    "        data = img.get_data()\n",
    "        tsnr = np.mean(data,axis=-1) / np.std(data,axis=-1)\n",
    "        tsnrs.append(tsnr)\n",
    "    mean_tsnr = np.mean(tsnrs,axis=0)\n",
    "    new_img = nb.Nifti1Image(mean_tsnr,affine=img.affine,header=img.header)\n",
    "    save_fn = os.path.join(ts_write_dir,'sub_rank_%s_tsnr_avg_over_tasks.nii'%(si))\n",
    "    nb.save(new_img,save_fn)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# create tsnr map for different participant groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "870759\n"
     ]
    }
   ],
   "source": [
    "sg={\n",
    "# 186: range(0,5),\n",
    "187: range(88,93),\n",
    "# 184: range(176,181)\n",
    "}\n",
    "\n",
    "for s in sg.keys():\n",
    "    tsnrs = []\n",
    "    for si in sg[s]:\n",
    "        fn = os.path.join(ts_write_dir,'sub_rank_%s_tsnr_avg_over_tasks.nii'%(si))\n",
    "        img = nb.load(fn)\n",
    "        tsnrs.append(img.get_data())\n",
    "    mean_tsnr = np.mean(tsnrs,axis=0)\n",
    "    print np.isnan(mean_tsnr).sum()\n",
    "    new_img = nb.Nifti1Image(mean_tsnr,affine=img.affine,header=img.header)\n",
    "    save_fn = os.path.join(ts_write_dir,'sub_rank_%d_tsnr_avg_over_tasks.nii'%(s))\n",
    "    nb.save(new_img,save_fn)\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# correlate tsnr with participant rank for each voxel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 91/91 [07:57<00:00,  5.18s/it]\n"
     ]
    }
   ],
   "source": [
    "tsnrs = []\n",
    "ranks = []\n",
    "for si in range(181):\n",
    "    fn = os.path.join(ts_write_dir,'sub_rank_%s_tsnr_avg_over_tasks.nii'%(si))\n",
    "    img = nb.load(fn)\n",
    "    tsnr = img.get_data()\n",
    "    tsnrs.append(tsnr)\n",
    "    ranks.append(np.ones_like(tsnr)*si)\n",
    "\n",
    "tsnrs = np.array(tsnrs)\n",
    "ranks = np.array(ranks)\n",
    "\n",
    "lenx = np.shape(tsnrs)[1]\n",
    "leny = np.shape(tsnrs)[2]\n",
    "lenz = np.shape(tsnrs)[3]\n",
    "\n",
    "corrs = np.zeros((lenx,leny,lenz))\n",
    "ps = np.zeros((lenx,leny,lenz))\n",
    "for x in tqdm(range(lenx)):\n",
    "    for y in range(leny):\n",
    "        for z in range(lenz):\n",
    "            r,p = stats.spearmanr(tsnrs[:,x,y,z],ranks[:,x,y,z])\n",
    "#             if p < .05:\n",
    "            corrs[x,y,z] = r\n",
    "#             else:\n",
    "#                 corrs[x,y,z] = np.nan\n",
    "            ps[x,y,z] = p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "new_img = nb.Nifti1Image(corrs,affine=img.affine,header=img.header)\n",
    "save_fn = os.path.join(ts_write_dir,'tsnr_rank_corrs.nii')\n",
    "nb.save(new_img,save_fn)\n",
    "\n",
    "new_img = nb.Nifti1Image(ps,affine=img.affine,header=img.header)\n",
    "save_fn = os.path.join(ts_write_dir,'tsnr_rank_corrs_p.nii')\n",
    "nb.save(new_img,save_fn)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [analysis]",
   "language": "python",
   "name": "Python [analysis]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}