A_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## import packages, define variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "warnings.simplefilter('ignore')\n",
    "import pandas as pd # for creation of data frames\n",
    "import numpy as np #\n",
    "import seaborn as sns\n",
    "import nibabel as nib # for loading niftis etc. as arrays\n",
    "import matplotlib.pyplot as plt # for plotting data\n",
    "#the output of plotting commands is displayed inline, directly below the code cell that produced it\n",
    "%matplotlib inline \n",
    "import os,glob\n",
    "from nilearn import plotting, input_data, image #for plotting & working with niftFMis\n",
    "\n",
    "import pathlib\n",
    "import pyxnat\n",
    "\n",
    "###################### update path!! #####################################\n",
    "base_path = '/root_dir/' ####\n",
    "##########################################################################\n",
    "\n",
    "#import qBOLD functions\n",
    "#import sys\n",
    "#sys.path.append(os.path.join(base_path + 'scripts/'))\n",
    "#import mqBOLD_functions as qB\n",
    "\n",
    "\n",
    "### define FSL directories ######################################################\n",
    "os.environ[\"FSLDIR\"]='/usr/share/fsl/5.0'\n",
    "os.environ[\"FSLOUTPUTTYPE\"]='NIFTI_GZ'\n",
    "os.environ[\"FSLTCLSH\"]='/usr/bin/tclsh'\n",
    "os.environ[\"FSLWISH\"]='/usr/bin/wish'\n",
    "os.environ[\"FSLMULTIFILEQUIT\"]=\"True\"\n",
    "os.environ[\"LD_LIBRARY_PATH\"]='/usr/share/fsl/5.0:/usr/lib/fsl/5.0'\n",
    "#################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = base_path + 'data/'\n",
    "derivatives_dir = data_dir + 'derivatives'\n",
    "raw_dir = data_dir + 'rawdata'\n",
    "results_dir = base_path + 'results'\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "\n",
    "## N=40\n",
    "sids = [19, 20, 21, 23, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 43, 44, 46, 47, 48, 49, 50, 51, 52, 54, 55, 58, 59, 60, 61, 63, 64, 65, 66, 67, 68]\n",
    "## N=30\n",
    "#sids_30 = [19, 20, 21, 23, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 43, 44, 46, 47, 48, 49, 50, 51, 52, 54, 55]\n",
    "## N=10\n",
    "#sids_10 = [58, 59, 60, 61, 63, 64, 65, 66, 67, 68]\n",
    "\n",
    "df_participants = pd.read_csv(data_dir + '/participants.tsv',sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate CMRO2 per subject and condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, ID in enumerate(sids): #loop over subjects\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "\n",
    "    #download & binarize brain mask, store in 'masker' object\n",
    "    BrMask_CSF = os.path.join(derivatives_dir, sub, 'anat',  sub + '_BrMsk_CSF.nii')\n",
    "    mask = image.math_img('img> 0.1', img = BrMask_CSF)\n",
    "    masker = input_data.NiftiMasker(mask_img=mask)\n",
    "    \n",
    "    # download O2 sat & hematocrit values from xnat\n",
    "    Hk = df_participants.query('participant_id==@sub')['Hct'].values\n",
    "    O2sat = df_participants.query('participant_id==@sub')['O2sat'].values\n",
    "        \n",
    "    if ID < 56: ## here, all 4 conditions were acquired        \n",
    "        conds = ['rest', 'control', 'mem', 'calc']\n",
    "        \n",
    "    if ID > 55:        \n",
    "        conds = ['control', 'calc'] \n",
    "        \n",
    "    for k, cond in enumerate(conds): \n",
    "\n",
    "        #  CBF and OEF niftis f\n",
    "        CBF = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-' + cond + '_space-T2_cbf.nii')\n",
    "        CBF_np = masker.fit_transform(CBF)\n",
    "        OEF = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-' + cond + '_space-T2_oef.nii')\n",
    "        OEF_np = masker.fit_transform(OEF)\n",
    "        \n",
    "        # upscale CBF by 25% because 4 background suppression pulses\n",
    "        CBF_upscaled = CBF_np / 0.75\n",
    "\n",
    "        #calculate CaO2 value per subject, O2sat as a fraction\n",
    "        #\n",
    "        CaO2 = 0.334 * float(Hk) * 55.6 * float(O2sat)/100;\n",
    "        \n",
    "        ## calculate CMRO2 parameter map\n",
    "        CMRO2 = CBF_upscaled * OEF_np * CaO2 /100;\n",
    "        \n",
    "        #store CMRO2 parameter map as nifti\n",
    "        CMRO2_vol = masker.inverse_transform(CMRO2)\n",
    "        ## here, CBV was only acquired in control\n",
    "        if ID < 56:\n",
    "            CMRO2_vol.to_filename(os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-' + cond + '_space-T2_desc-orig_cmro2.nii'))\n",
    "        ## here, CBV correction was applied: CBV was also acquired in calc\n",
    "        if ID > 55 and cond == 'calc':\n",
    "            CMRO2_vol.to_filename(os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-' + cond + '_space-T2_desc-CBV_cmro2.nii'))        \n",
    "        if ID > 55 and cond == 'control':\n",
    "            CMRO2_vol.to_filename(os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-' + cond + '_space-T2_desc-orig_cmro2.nii')) \n",
    "            \n",
    "        ## plot image\n",
    "        #plotting.plot_img(CMRO2_vol, cmap='jet', colorbar=True, vmin=10, vmax=300, title=sub + ': ' + cond)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CBV calc-correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### create 'new' CBV calc images ###\n",
    "\n",
    "df_Glasser = pd.read_csv(derivatives_dir + '/CBV_correction_factor.tsv',sep='\\t')\n",
    "correction_factor_percchange = df_Glasser['CBV_percchange']\n",
    "\n",
    "for i,ID in enumerate(sids_30):\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "\n",
    "    ## download CBV control map\n",
    "    CBV_control = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-control_space-T2_cbv.nii')\n",
    "    CBV_control_img = nib.load(CBV_control)\n",
    "    CBV_control_arr = np.array(CBV_control_img.dataobj)\n",
    "    \n",
    "    ## download Glasser360 parcellation\n",
    "    Glasser_T2space = os.path.join(derivatives_dir, 'Glasser360_T2space.nii.gz')\n",
    "    if ID > 23:\n",
    "        Glasser_T2space = os.path.join(derivatives_dir, 'Glasser360_T2space_35vol.nii.gz')\n",
    "    Glasser_img=nib.load(Glasser_T2space)\n",
    "    Glasser_arr=np.array(Glasser_img.dataobj)\n",
    "       \n",
    "    ## CBV control img is the baseline, replace voxels that are in Glasser parcels with adjusted CBV values\n",
    "    CBV_calc_corrected = CBV_control_arr \n",
    "    for r, region in enumerate(range(1, 361)): \n",
    "        correction_factor = correction_factor_percchange/100 ##from percent to ratio\n",
    "        CBV_calc_corrected[Glasser_arr==region] = CBV_control_arr[Glasser_arr==region] + (correction_factor[r]*CBV_control_arr[Glasser_arr==region])\n",
    "    \n",
    "\n",
    "    ## turn into nifti, display\n",
    "    CBV_calc_corrected_img = nib.Nifti1Image(CBV_calc_corrected, CBV_control_img.affine)\n",
    "    T1w_T2 = os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_T1w.nii')\n",
    "\n",
    "    #plotting.plot_img(CBV_calc_corrected_img,bg_img=T1w_T2, cmap='jet', colorbar=True, threshold=0.1, vmax=10, black_bg=True, \n",
    "    #                             title = 'CBV calc corrected for sub ' + subj_id)\n",
    "    #plotting.plot_img(CBV_calc_corrected_img, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=10, threshold=0.1,\n",
    "    #               title = 'CBV calc corrected for sub ' + subj_id)\n",
    "    #plotting.plot_img(CBV_control_img, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=10, threshold=0.1,\n",
    "    #               title = 'CBV control for sub ' + subj_id)\n",
    "    \n",
    "    ## save \n",
    "    nib.save(CBV_calc_corrected_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_cbv.nii'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### recalculate OEF and CMRO2 for calc ###\n",
    "\n",
    "for i,ID in enumerate(sids_30):\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "    \n",
    "    # download O2 sat & hematocrit values from xnat\n",
    "    Hct = df_participants.query('participant_id==@sub')['Hct'].values\n",
    "    O2sat = df_participants.query('participant_id==@sub')['O2sat'].values\n",
    "        \n",
    "    ##############################\n",
    "    ## calculate adjusted rOEF ##\n",
    "    #############################\n",
    "    \n",
    "    \n",
    "    R2prime_calc =  os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_R2prime.nii')\n",
    "    R2prime_calc_nii = nib.load(R2prime_calc)\n",
    "    R2prime_calc_arr = np.array(R2prime_calc_nii.dataobj)\n",
    "    \n",
    "    ## co-register CBV and T2, do not have the same number of slices for ID<26! \n",
    "    CBV_calc_corrected =  os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_cbv.nii')\n",
    "    CBV_calc_corrected_nii = nib.load(CBV_calc_corrected)\n",
    "    CBV_calc_corrected_arr = np.array(CBV_calc_corrected_nii.dataobj)\n",
    "    \n",
    "    BrMask_CSF = os.path.join(derivatives_dir, sub, 'anat',  sub + '_BrMsk_CSF.nii')\n",
    "    BrMsk_CSF_nii = nib.load(BrMsk_CSF)\n",
    "    BrMsk_CSF_arr = np.array(BrMsk_CSF_nii.dataobj)\n",
    "      \n",
    "    ## here, the brain mask has too many slices\n",
    "    T1w_T2 = os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_T1w.nii')\n",
    "    T1w_T2_brain = os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_desc-brain_T1w.nii.gz')\n",
    "    ! bet {T1w_T2} {T1w_T2_brain}\n",
    "    BrMask_CSF_sliced = os.path.join(derivatives_dir, sub, 'anat',  sub + '_BrMsk_CSF_resliced.nii.gz')\n",
    "    ! flirt -in {BrMask_CSF} -ref {T1w_T2_brain} -out {BrMask_CSF_sliced} -dof 6 \n",
    "    plotting.plot_img(BrMask_CSF_sliced, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=350, vmin = 10, threshold=0.1,\n",
    "                   title = 'CMRO2 calc corrected: ' + sub)\n",
    "    \n",
    "    BrMsk_CSF_nii = nib.load(BrMask_CSF_sliced)\n",
    "    BrMsk_CSF_arr = np.array(BrMsk_CSF_nii.dataobj)\n",
    "\n",
    "        \n",
    "    mask = BrMsk_CSF_arr>0.5;\n",
    "    \n",
    "    # [Hz], physiologic constant derived from literature\n",
    "    # 0.85: factor for small-vessel hematocrit\n",
    "    C = 4/3*267.61918*np.pi*0.264*(Hct/100)*0.85*3; # 317 Hz\n",
    "\n",
    "    rCBV = CBV_calc_corrected_arr/100.0; #CBV fraction, e.g. CBV_WM = 0.015\n",
    "    rCBV[rCBV < 0] = 0;\n",
    "    eps=2.2204e-16 #to not divide by zero\n",
    "    rOEF = (R2prime_calc_arr / (C*rCBV+eps)) * mask\n",
    "\n",
    "    # reasonable boundary to avoid unreasonably high hot spots\n",
    "    rOEFmax_DSC = 5*np.nanmedian(rOEF[mask]);\n",
    "    if rOEFmax_DSC < 1.5:\n",
    "            rOEFmax_DSC = 1.5 # FIXED Threshold\n",
    "    rOEF[rOEF > rOEFmax_DSC] = rOEFmax_DSC\n",
    "    rOEF = rOEF *mask\n",
    "\n",
    "    rOEF_calc_corrected_img = nib.Nifti1Image(rOEF, CBV_calc_corrected_nii.affine)\n",
    "    nib.save(rOEF_calc_corrected_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_desc-CBV_oef.nii'))\n",
    "\n",
    "\n",
    "    ##############################\n",
    "    ## calculate adjusted CMRO2 ##\n",
    "    ##############################\n",
    "    \n",
    "    \n",
    "    # download CSF and OEF niftis from xnat\n",
    "    CBF = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_cbf.nii')\n",
    "    CBF_nii = nib.load(CBF)\n",
    "    CBF_np= np.array(CBF_nii.dataobj)\n",
    "    OEF_corrected = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_desc-CBV_oef.nii')\n",
    "    OEF_nii = nib.load(OEF_corrected)\n",
    "    OEF_np= np.array(OEF_nii.dataobj)\n",
    "    \n",
    "    # upscale CSF by 25% because of scanner\n",
    "    CBF_upscaled = CBF_np / 0.75 # upscale CSF by 25% because of scanner\n",
    "\n",
    "    #calculate CaO2 value per subject, O2sat as a fraction\n",
    "    CaO2 = 0.334 * Hct * 55.6 * float(O2sat)/100;\n",
    "        \n",
    "    ## calculate CMRO2 parameter map\n",
    "    CMRO2 = CBF_upscaled * OEF_np * CaO2 /100;\n",
    "    CMRO2 = CMRO2 * mask\n",
    "        \n",
    "    #store CMRO2 parameter map as nifti\n",
    "    CMRO2_calc_corrected_img = nib.Nifti1Image(CMRO2, CBV_calc_corrected_nii.affine)\n",
    "\n",
    "    nib.save(CMRO2_calc_corrected_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_desc-CBV_cmro2.nii.gz'))\n",
    "\n",
    "    ## plot output\n",
    "    T1w_T2 = os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_T1w.nii')\n",
    "    plotting.plot_img(CMRO2_calc_corrected_img, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=400, vmin = 10, threshold=0.1,\n",
    "                   title = 'CMRO2 calc corrected: ' + sub)\n",
    "    CMRO2_calc_img = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_desc-orig_cmro2.nii')\n",
    "    plotting.plot_img(CMRO2_calc_img, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=400, vmin = 10, threshold=0.1,\n",
    "                   title = 'CMRO2 orig: ' + sub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### store BOLD percchange niftis (desc-preproc instead of filtered_func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## coregister func to T2 map\n",
    "for i, ID in enumerate(sids): #loop over subjects\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "    sub_ID = \"p{:03d}\".format(ID) \n",
    "\n",
    "    shutil.copy('/home/tumnic/sepp/neuroenergeticslab_notebooks/erc-wp1/tmp/'+sub_ID+'/fMRI/fmriprep/sub-'+sub_ID+'/func/sub-'+sub_ID+'_task-locmem_space-T1w_desc-preproc_bold.nii.gz', os.path.join(derivatives_dir, sub, 'func',  sub + '_task-all_space-T1w_desc-preproc_bold.nii.gz'))\n",
    "    \n",
    "    T1w_T1space = os.path.join(derivatives_dir, sub, 'anat',sub + '_desc-fmriprep_T1w.nii.gz')\n",
    "    T1w_T2space = os.path.join(derivatives_dir, sub, 'qmri',sub + '_space-T2_T1w.nii')\n",
    "    T1w_T2space_brain = os.path.join(derivatives_dir, sub, 'qmri',sub + '_space-T2_T1w_brain.nii.gz')\n",
    "    T1w_to_T2space = os.path.join(derivatives_dir, sub, 'qmri',sub + '_T1w_to_T2space.mat')\n",
    "    \n",
    "    func_T1space = os.path.join(derivatives_dir, sub, 'func',  sub + '_task-all_space-T1w_desc-preproc_bold.nii.gz')    \n",
    "    func_T2space = os.path.join(derivatives_dir, sub, 'func',  sub + '_task-all_space-T2_desc-preproc_bold.nii.gz')\n",
    "    \n",
    "    ## brain-extract T1w T2 space\n",
    "    ! bet {T1w_T2space} {T1w_T2space_brain} -f 0.3\n",
    "    plotting.plot_img(T1w_T2space_brain,  bg_img =T1w_T2space, cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), threshold = 1, black_bg=True, title = 'T1w_brain : ' + sub)\n",
    "    \n",
    "    ## coregister T1 to T2 space \n",
    "    ! epi_reg --epi={func_T1space} --t1={T1w_T2space} --t1brain={T1w_T2space_brain} --out={func_T2space}\n",
    "\n",
    "    \n",
    "    func_T2space_slice = image.index_img(func_T2space, 0) ## first volume\n",
    "    plotting.plot_img(func_T2space_slice,  bg_img =T1w_T2space, cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), threshold = 10000, black_bg=True, title = 'func : ' + sub)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, ID in enumerate(sids): #loop over subjects\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "    \n",
    "    if ID < 56:             \n",
    "        contrasts = ['calccontrol', 'memcontrol', 'calcrest']\n",
    "        #contrasts = ['memcontrol']\n",
    "                    \n",
    "        func_T2space = os.path.join(derivatives_dir, sub, 'func',  sub + '_task-all_space-T2_desc-preproc_bold.nii.gz')\n",
    "        func_img = nib.load(func_T2space)\n",
    "        func_arr = np.array(func_img.dataobj) \n",
    "\n",
    "        ##taking HRF into account, i.e. only take time-points after 5TR = 7sec\n",
    "        rest = func_arr[:, :, :,list(range(5, 25)) + list(range(130, 150)) + list(range(230, 250)) + list(range(330, 350))] #rest\n",
    "        baseline_rest = np.nanmedian(rest, axis=3) #median across time-points\n",
    "\n",
    "        control = func_arr[:, :, :, list(range(80, 100)) + list(range(205, 225)) + list(range(280, 300)) + list(range(355, 375))]\n",
    "        baseline_control = np.nanmedian(control, axis=3) #median across time-points\n",
    "\n",
    "                \n",
    "        ## get %change per time-point\n",
    "        percchange_control=np.zeros((func_arr.shape))\n",
    "        percchange_rest=np.zeros((func_arr.shape))\n",
    "        percchange_mem = np.zeros((func_arr.shape))\n",
    "        for i in range(0, func_arr.shape[3]):\n",
    "            percchange_control[:, :, :, i] = 100*(func_arr[:, :, :, i] - baseline_control) / baseline_control\n",
    "            percchange_rest[:, :, :, i] = 100*(func_arr[:, :, :, i] - baseline_rest) / baseline_rest\n",
    "\n",
    "        for l, con in enumerate(contrasts): \n",
    "\n",
    "            if con == 'calccontrol':\n",
    "                percchange = np.nanmedian(percchange_control[:, :, :,list(range(55, 75)) + list(range(155, 175)) + list(range(255, 275)) + list(range(305, 325))], axis=3)\n",
    "\n",
    "            if con == 'memcontrol':\n",
    "                percchange = np.nanmedian(percchange_control[:, :, :,list(range(30, 50)) + list(range(105, 125)) + list(range(180, 200)) + list(range(380, 400))], axis=3)\n",
    "\n",
    "            if con == 'calcrest':\n",
    "                percchange = np.nanmedian(percchange_rest[:, :, :,list(range(55, 75)) + list(range(155, 175)) + list(range(255, 275)) + list(range(305, 325))], axis=3)\n",
    "\n",
    "            ## save nifti\n",
    "            T1w_T2 = nib.load(os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_T1w.nii'))\n",
    "            percchange_img = nib.Nifti1Image(percchange, T1w_T2.affine)\n",
    "            plotting.plot_img(percchange_img,  cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), colorbar=True, black_bg=True, threshold=0, vmin=0, vmax=5, title = 'percchange : ' + sub)\n",
    "\n",
    "            nib.save(percchange_img, os.path.join(derivatives_dir, sub, 'func',  sub +  '_task-'+con+'_space-T2_BOLD_percchange.nii.gz'))\n",
    "\n",
    "\n",
    "    if ID > 55:\n",
    "        con = 'calccontrol'\n",
    "        \n",
    "        func_T2space = os.path.join(derivatives_dir, sub, 'func',  sub +  '_task-all_space-T2_desc-preproc_bold.nii.gz')\n",
    "        func_img = nib.load(func_T2space)\n",
    "        func_arr = np.array(func_img.dataobj)             \n",
    "\n",
    "        ##taking HRF into account, i.e. only take time-points after 5TR = 7sec\n",
    "        control = func_arr[:, :, :, list(range(5, 25)) + list(range(55, 75)) + list(range(105, 125)) + list(range(155, 175))]\n",
    "        baseline_control = np.nanmedian(control, axis=3) #median across time-points\n",
    "\n",
    "        ## get %change per time-point\n",
    "        percchange_control=np.zeros((func_arr.shape))\n",
    "        for i in range(0, func_arr.shape[3]):\n",
    "            percchange_control[:, :, :, i] = 100*(func_arr[:, :, :, i] - baseline_control) / baseline_control\n",
    "\n",
    "        percchange = np.nanmedian(percchange_control[:, :, :,list(range(30, 50)) + list(range(80, 100)) + list(range(130, 150)) + list(range(180, 200))], axis=3)\n",
    "\n",
    "        ## save nifti\n",
    "        T1w_T2 = nib.load(os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_T1w.nii'))\n",
    "        percchange_img = nib.Nifti1Image(percchange, T1w_T2.affine)\n",
    "        nib.save(percchange_img, os.path.join(derivatives_dir, sub, 'func',  sub +  '_task-'+con+'_space-T2_BOLD_percchange.nii.gz'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BOLD percchange MNI space based on fmriprep output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, ID in enumerate(sids): #loop over subjects\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "    \n",
    "    if ID < 56:             \n",
    "        contrasts = ['calccontrol', 'memcontrol', 'calcrest']\n",
    "        #contrasts = ['calcrest']\n",
    "                    \n",
    "        func_T2space = os.path.join(derivatives_dir, sub, 'func',  sub +  'task-all_space-MNI152_res-2_desc-preproc_bold.nii.gz')\n",
    "        func_img = nib.load(func_T2space)\n",
    "        func_arr = np.array(func_img.dataobj) \n",
    "\n",
    "        ##taking HRF into account, i.e. only take time-points after 5TR = 7sec\n",
    "        rest = func_arr[:, :, :,list(range(5, 25)) + list(range(130, 150)) + list(range(230, 250)) + list(range(330, 350))] #rest\n",
    "        baseline_rest = np.nanmedian(rest, axis=3) #median across time-points\n",
    "\n",
    "        control = func_arr[:, :, :, list(range(80, 100)) + list(range(205, 225)) + list(range(280, 300)) + list(range(355, 375))]\n",
    "        baseline_control = np.nanmedian(control, axis=3) #median across time-points\n",
    "\n",
    "                \n",
    "        ## get %change per time-point\n",
    "        percchange_control=np.zeros((func_arr.shape))\n",
    "        percchange_rest=np.zeros((func_arr.shape))\n",
    "        percchange_mem = np.zeros((func_arr.shape))\n",
    "        for i in range(0, func_arr.shape[3]):\n",
    "            percchange_control[:, :, :, i] = 100*(func_arr[:, :, :, i] - baseline_control) / baseline_control\n",
    "            percchange_rest[:, :, :, i] = 100*(func_arr[:, :, :, i] - baseline_rest) / baseline_rest\n",
    "\n",
    "        for l, con in enumerate(contrasts): \n",
    "\n",
    "            if con == 'calccontrol':\n",
    "                percchange = np.nanmedian(percchange_control[:, :, :,list(range(55, 75)) + list(range(155, 175)) + list(range(255, 275)) + list(range(305, 325))], axis=3)\n",
    "\n",
    "            if con == 'memcontrol':\n",
    "                percchange = np.nanmedian(percchange_control[:, :, :,list(range(30, 50)) + list(range(105, 125)) + list(range(180, 200)) + list(range(380, 400))], axis=3)\n",
    "\n",
    "            if con == 'calcrest':\n",
    "                percchange = np.nanmedian(percchange_rest[:, :, :,list(range(55, 75)) + list(range(155, 175)) + list(range(255, 275)) + list(range(305, 325))], axis=3)\n",
    "\n",
    "            ## save nifti\n",
    "            T1w_T2 = nib.load(os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz'))\n",
    "            percchange_img = nib.Nifti1Image(percchange, T1w_T2.affine)\n",
    "            #plotting.plot_img(percchange_img,  cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), colorbar=True, black_bg=True, threshold=0, vmin=0, vmax=5, title = 'percchange : ' + sub)\n",
    "\n",
    "            nib.save(percchange_img, os.path.join(derivatives_dir, sub, 'func',  sub +  '_task-'+con+'_space-MNI152_desc-fmriprep_BOLD_percchange.nii.gz'))\n",
    "\n",
    "\n",
    "    if ID > 55:\n",
    "        con = 'calccontrol'\n",
    "        \n",
    "        func_T2space = os.path.join(derivatives_dir, sub, 'func',  sub +  'task-all_space-MNI152_res-2_desc-preproc_bold.nii.gz')\n",
    "        func_img = nib.load(func_T2space)\n",
    "        func_arr = np.array(func_img.dataobj)             \n",
    "\n",
    "        ##taking HRF into account, i.e. only take time-points after 5TR = 7sec\n",
    "        control = func_arr[:, :, :, list(range(5, 25)) + list(range(55, 75)) + list(range(105, 125)) + list(range(155, 175))]\n",
    "        baseline_control = np.nanmedian(control, axis=3) #median across time-points\n",
    "\n",
    "        ## get %change per time-point\n",
    "        percchange_control=np.zeros((func_arr.shape))\n",
    "        for i in range(0, func_arr.shape[3]):\n",
    "            percchange_control[:, :, :, i] = 100*(func_arr[:, :, :, i] - baseline_control) / baseline_control\n",
    "\n",
    "        percchange = np.nanmedian(percchange_control[:, :, :,list(range(30, 50)) + list(range(80, 100)) + list(range(130, 150)) + list(range(180, 200))], axis=3)\n",
    "\n",
    "        ## save nifti\n",
    "        T1w_T2 = nib.load(os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz'))\n",
    "        percchange_img = nib.Nifti1Image(percchange, T1w_T2.affine)\n",
    "        nib.save(percchange_img, os.path.join(derivatives_dir, sub, 'func',  sub +  '_task-'+con+'_space-MNI152_desc-fmriprep_BOLD_percchange.nii.gz'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### semi-quantitative CMRO2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys \n",
    "#!{sys.prefix}/bin/pip install -q nbformat \n",
    "import nbformat "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## run qBOLD toolbox\n",
    "%run qBOLD_fun.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## which contrast/baseline do we want to look at? ##\n",
    "#############################################################\n",
    "baseline='rest'\n",
    "#baseline='control'\n",
    "#conds=['calc']\n",
    "#conds=['mem']\n",
    "\n",
    "### calculate based on CBV-corrected calc or without CBV-correction?\n",
    "mode = 'corrected'\n",
    "#mode = 'noCBV'\n",
    "#############################################################\n",
    "\n",
    "coords=(-30, -15, 0, 15, 30, 45, 60, 75, 90)\n",
    "\n",
    "TE=30\n",
    "\n",
    "for i, ID in enumerate(sids): #loop over subjects\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "    \n",
    "    # download O2 sat & hematocrit values from xnat\n",
    "    Hct = df_participants.query('participant_id==@sub')['Hct'].values\n",
    "    O2sat = df_participants.query('participant_id==@sub')['O2sat'].values\n",
    "    \n",
    "    ## load R2` and %BOLD change\n",
    "    R2s_baseline =  os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+baseline+'_space-T2_R2prime.nii')\n",
    "    R2s_baseline_nii = nib.load(R2s_baseline)\n",
    "    R2s_baseline_arr = np.array(R2s_baseline_nii.dataobj)\n",
    "   # plotting.plot_img(R2s_baseline, cmap='jet', cut_coords=coords, display_mode='z',colorbar = True, vmin = -20, vmax=20, title=subj_id + ' R2s baseline')\n",
    "    \n",
    "    if ID < 56:\n",
    "        conds=['mem', 'calc']\n",
    "        #conds=['calc']\n",
    "    if ID > 55:\n",
    "        conds=['calc']\n",
    "        \n",
    "    for cond in conds:\n",
    "        con = cond + baseline\n",
    "        BOLD = os.path.join(derivatives_dir, sub, 'func',  sub +  '_task-'+con+'_space-T2_BOLD_percchange.nii.gz') ## input: desc-preproc_bold.nii.gz from fmriprep\n",
    "    \n",
    "        #plotting.plot_img(BOLD, cmap='jet', cut_coords=coords, display_mode='z',colorbar = True, vmin = -2, vmax=2, title=subj_id + ' BOLD percchange')\n",
    "        BOLD_nii = nib.load(BOLD)\n",
    "        BOLD_arr = np.array(BOLD_nii.dataobj)\n",
    "\n",
    "        ########################\n",
    "        ## calculate R2'     ###\n",
    "        ########################\n",
    "\n",
    "        ##recalcalute task R2' based on R2' baseline and BOLD change\n",
    "        R2s_Fujita = FujitaR2ptask_nii(R2s_baseline_arr, BOLD_arr, TE)\n",
    "        ## turn into nifti and save\n",
    "        R2s_Fujita_nii = nib.Nifti1Image(R2s_Fujita, R2s_baseline_nii.affine)\n",
    "        nib.save(R2s_Fujita_nii, os.path.join(derivatives_dir, sub, 'qmri',   sub + '_task-'+cond+'_space-T2_desc-semi-quant_R2prime.nii.gz'))\n",
    "        #plotting.plot_img(R2s_Fujita_nii, cmap='jet', cut_coords=coords, display_mode='z',colorbar = True, vmin = -20, vmax=20, title=subj_id + ' R2s Fujita')\n",
    "\n",
    "        ###################\n",
    "        ## calculate OEF ##\n",
    "        ###################\n",
    "        ## download corrected CBV map\n",
    "        if cond == 'calc': \n",
    "            if mode == 'corrected':## for calc, we take the corrected CBV values!\n",
    "                CBV =  os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-calc_space-T2_cbv.nii')\n",
    "            if mode == 'noCBV': ## if not corrected, we take baseline CBV values\n",
    "                CBV =  os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-control_space-T2_cbv.nii')\n",
    "        if cond == 'mem': # for the other conditions, we take the CBV from the control task\n",
    "            CBV =  os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-control_space-T2_cbv.nii')\n",
    "        CBV_nii = nib.load(CBV)\n",
    "        CBV_arr = np.array(CBV_nii.dataobj)\n",
    "\n",
    "        ##download brain mask\n",
    "        BrMsk_CSF = os.path.join(derivatives_dir, sub, 'anat',  sub + '_BrMsk_CSF.nii')\n",
    "        BrMsk_CSF_nii = nib.load(BrMsk_CSF)\n",
    "        BrMsk_CSF_arr = np.array(BrMsk_CSF_nii.dataobj)\n",
    "        #plotting.plot_img(BrMsk_CSF_nii, cmap='jet', threshold=0.5, cut_coords=coords, display_mode='z',colorbar = True, vmin = -1, vmax=1, title=subj_id + ' mask')\n",
    "\n",
    "        mask = BrMsk_CSF_arr>0.5;\n",
    "\n",
    "        # [Hz], physiologic constant derived from literature\n",
    "        # 0.85: factor for small-vessel hematocrit\n",
    "        C = 4/3*267.61918*np.pi*0.264*(Hct/100)*0.85*3; # 317 Hz\n",
    "\n",
    "        rCBV = CBV_arr/100.0; #CBV fraction, e.g. CBV_WM = 0.015\n",
    "        rCBV[rCBV < 0] = 0;\n",
    "        eps=2.2204e-16 #to not divide by zero\n",
    "        rOEF = (R2s_Fujita / (C*rCBV+eps)) * mask\n",
    "\n",
    "        # reasonable boundary to avoid unreasonably high hot spots\n",
    "        #rOEFmax_DSC = 5*np.nanmedian(rOEF[mask]);\n",
    "        rOEFmax_DSC = 5*np.nanmedian(rOEF[mask]);\n",
    "        if rOEFmax_DSC < 1.5:\n",
    "                rOEFmax_DSC = 1.5 # FIXED Threshold\n",
    "        rOEF[abs(rOEF) > rOEFmax_DSC] = rOEFmax_DSC\n",
    "        rOEF = rOEF *mask\n",
    "\n",
    "        rOEF_img = nib.Nifti1Image(rOEF, R2s_baseline_nii.affine)\n",
    "        if (mode == 'corrected') and (cond =='calc'):\n",
    "            nib.save(rOEF_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant-corrected_oef.nii'))\n",
    "        if (mode == 'noCBV') and (cond =='calc'):\n",
    "            nib.save(rOEF_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_oef.nii'))\n",
    "        if (mode == 'noCBV') or (cond !='calc'):\n",
    "            nib.save(rOEF_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_oef.nii'))\n",
    "\n",
    "        #####################\n",
    "        ## calculate CMRO2 ##\n",
    "        #####################\n",
    "\n",
    "        # download CBF and OEF niftis from xnat\n",
    "        CBF = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_space-T2_cbf.nii')\n",
    "        #plotting.plot_img(CBF, cmap='jet', cut_coords=coords, display_mode='z',colorbar = True, vmin = -50, vmax=50, title=subj_id + ' CBF')\n",
    "\n",
    "        CBF_nii = nib.load(CBF)\n",
    "        CBF_np= np.array(CBF_nii.dataobj)\n",
    "        # upscale CSF by 25% because of scanner\n",
    "        CBF_upscaled = CBF_np / 0.75 # upscale CSF by 25% because of scanner\n",
    "        \n",
    "        if mode == 'corrected' and cond == 'calc':\n",
    "            OEF_Fujita = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant-corrected_oef.nii')\n",
    "        if mode == 'noCBV' and cond =='calc':\n",
    "            OEF_Fujita = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_oef.nii')\n",
    "        if mode == 'noCBV' or cond != 'calc':\n",
    "            OEF_Fujita = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_oef.nii')\n",
    "\n",
    "        OEF_nii = nib.load(OEF_Fujita)\n",
    "        OEF_np= np.array(OEF_nii.dataobj)\n",
    "\n",
    "        #calculate CaO2 value per subject, O2sat as a fraction\n",
    "        CaO2 = 0.334 * Hct * 55.6 * float(O2sat)/100;\n",
    "\n",
    "        ## calculate CMRO2 parameter map\n",
    "        CMRO2 = CBF_upscaled * OEF_np * CaO2 /100;\n",
    "        CMRO2 = CMRO2 * mask\n",
    "\n",
    "        #store CMRO2 parameter map as nifti\n",
    "        CMRO2_img = nib.Nifti1Image(CMRO2, CBV_nii.affine)\n",
    "        if mode == 'corrected' and cond == 'calc':\n",
    "            #nib.save(CMRO2_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant-corrected_cmro2.nii.gz'))\n",
    "            nib.save(CMRO2_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant-corrected_cmro2.nii.gz'))            \n",
    "        if mode == 'noCBV' and cond =='calc':\n",
    "            #nib.save(CMRO2_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_cmro2.nii.gz'))                     \n",
    "            nib.save(CMRO2_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_cmro2.nii.gz'))                     \n",
    "        if mode == 'noCBV' or cond != 'calc':\n",
    "            #nib.save(CMRO2_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_cmro2.nii.gz'))                     \n",
    "            nib.save(CMRO2_img, os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_base-'+baseline+'_space-T2_desc-semi-quant_cmro2.nii.gz'))                     \n",
    "            \n",
    "        ## plot output\n",
    "        T1w_T2 = os.path.join(derivatives_dir, sub, 'anat',  sub + '_space-T2_T1w.nii')\n",
    "        plotting.plot_img(CMRO2_img, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=400, vmin = 10, threshold=0.1,\n",
    "                   cut_coords=coords, title = 'CMRO2 semi-quant: '+ cond + ' ' + sub)\n",
    "        if ID < 56:\n",
    "            CMRO2_calc_img = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_space-T2_desc-orig_cmro2.nii') ## this is the original cmro2 image\n",
    "        if ID > 55:\n",
    "            CMRO2_calc_img = os.path.join(derivatives_dir, sub, 'qmri',  sub + '_task-'+cond+'_space-T2_desc-CBV_cmro2.nii') ## this is the original cmro2 image\n",
    "            \n",
    "        plotting.plot_img(CMRO2_calc_img, bg_img=T1w_T2, cmap='jet',display_mode='z', colorbar=True, black_bg=True, vmax=400, vmin = 10, threshold=0.1,\n",
    "                   cut_coords=coords, title = 'CMRO2 orig: ' + cond + ' ' + sub)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### normalize parameter maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    ". /etc/fsl/5.0/fsl.sh\n",
    ". /etc/ants/ants.sh\n",
    "\n",
    "parameters=\"cmro2 R2prime T2Smap oef cbv\"\n",
    "parameters=\"cmro2 cbf oef\"\n",
    "\n",
    "sids=\"19 20 21 23 26 27 28 29 30 31 32 34 33 35 36 37 38 39 40 43 44 46 47 48 49 50 51 52 54 55 58 59 60 61 63 64 65 66 67 68\"\n",
    "#sids=\"35 36 37 38 39 40 43 44 46 47 48 49 50 51 52 54 55 58 59 60 61 63 64 65 66 67 68\"\n",
    "\n",
    "#sids=\"19 20\"\n",
    "\n",
    "for ID in $sids; \n",
    "do\n",
    "    ##define directories\n",
    "    SID=p0${ID}\n",
    "    echo $SID\n",
    "    base_path=\"/home/tumnic/sepp/neuroenergeticslab_notebooks/erc-wp1/manuscript\" \n",
    "    locfolder=\"${base_path}/data/derivatives/sub-${SID}\"\n",
    "    \n",
    "    ##define input\n",
    "    ## output from fmriprep\n",
    "    T1w=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T1w.nii.gz\"\n",
    "    br_mask=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_brain_mask.nii.gz\"\n",
    "    transform_to_standard=\"${locfolder}/anat/sub-${SID}_from-T1w_to-MNI152NLin6Asym_mode-image_xfm.h5\"\n",
    "    #echo $transform_to_standard\n",
    "    T1w_MNIspace=\"${locfolder}/anat/sub-${SID}_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz\"\n",
    "    #echo $T1w_MNIspace\n",
    "    ## output from qmri\n",
    "    T1w_T2space=\"${locfolder}/qmri/sub-${SID}_space-T2_T1w.nii\"\n",
    "    \n",
    "    ##define output\n",
    "    T1w_brain=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T1w_brain.nii.gz\"\n",
    "    T2_to_T1w=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T2_to_T1w.mat\"\n",
    "    \n",
    "    ##brain-extract T1w\n",
    "    echo \"apply brain mask for subject ${SID} (BET)\"\n",
    "    #fslmaths $T1w -mas $br_mask $T1w_brain\n",
    "    #bet $T1w $T1w_brain\n",
    "    \n",
    "    ##get T2 to T1w transformation mat\n",
    "    #flirt -in $T1w_T2space -ref $T1w_brain -omat $T2_to_T1w -dof 6\n",
    "\n",
    "    \n",
    "    for par in $parameters; do    \n",
    "        if [ \"$ID\" -lt \"56\" ]; then\n",
    "        if [ par != \"cbv\" ]; then\n",
    "            conditions=\"rest control mem calc\"\n",
    "            #conditions=\"mem\"\n",
    "        fi\n",
    "        fi\n",
    "        if [ \"$ID\" -gt \"55\" ]; then\n",
    "            conditions=\"control calc\"\n",
    "        fi\n",
    "        if [ $par == \"cbv\" ]; then\n",
    "        if [ \"$ID\" -lt \"56\" ]; then\n",
    "            conditions=\"control\"\n",
    "        fi\n",
    "        fi\n",
    "      \n",
    "        \n",
    "        for cond in $conditions; do\n",
    "            echo $cond\n",
    "        \n",
    "            if [ $par != \"cmro2\" ]; then\n",
    "                parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-T2_${par}.nii\"\n",
    "                ##define outputs\n",
    "                out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-T1w_${par}.nii.gz\"\n",
    "                out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-MNI152_${par}.nii.gz\"\n",
    "            fi\n",
    "            \n",
    "            if [ $par == \"cmro2\" ] && [ $cond != \"calc\" ]; then\n",
    "                parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-T2_desc-orig_${par}.nii\"\n",
    "                out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${cond}_T1w-space_desc-orig_${par}.nii.gz\"\n",
    "                out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-MNI152_desc-orig_${par}.nii.gz\"\n",
    "            fi\n",
    "            if [ $par == \"cmro2\" ] && [ $cond == \"calc\" ] && [ $ID < 56 ]; then ##non-CBV corrected, but only for ID<56\n",
    "                parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-T2_desc-orig_${par}.nii\"\n",
    "                out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${cond}_T1w-space_desc-orig_${par}.nii.gz\"\n",
    "                out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-MNI152_desc-orig_${par}.nii.gz\"\n",
    "            fi\n",
    "            \n",
    "            ##apply transformation from lowres T2 to highres T1 space for parameter maps\n",
    "            flirt -in $parameter_map -ref $T1w_brain -applyxfm -init $T2_to_T1w -out $out_T1w\n",
    "\n",
    "            ##apply transformation (saved by fmriprep) from T1 to standard space\n",
    "            antsApplyTransforms -d 3 --float 1 --verbose 1 -i $out_T1w -r $T1w_MNIspace -t $transform_to_standard -o $out_MNI\n",
    "            \n",
    "            if [ $par == \"cmro2\" ] && [ $cond == \"calc\" ]; then ## all CBV-corrected\n",
    "                parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-T2_desc-CBV_${par}.nii\" \n",
    "                out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${cond}_T1w-space_desc-CBV_${par}.nii.gz\"\n",
    "                out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${cond}_space-MNI152_desc-CBV_${par}.nii.gz\"\n",
    "            \n",
    "                ##apply transformation from lowres T2 to highres T1 space for parameter maps\n",
    "                flirt -in $parameter_map -ref $T1w_brain -applyxfm -init $T2_to_T1w -out $out_T1w\n",
    "\n",
    "                ##apply transformation (saved by fmriprep) from T1 to standard space\n",
    "                antsApplyTransforms -d 3 --float 1 --verbose 1 -i $out_T1w -r $T1w_MNIspace -t $transform_to_standard -o $out_MNI\n",
    "            fi\n",
    "\n",
    "            \n",
    "        done\n",
    "    done\n",
    "done\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## show normalized parameter maps\n",
    "parameters = ['cbf', 'cmro2', 'oef']\n",
    "\n",
    "for par in parameters:\n",
    "    for ID in sids:\n",
    "        sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "        print(sub)\n",
    "        sub_ID = \"p{:03d}\".format(ID) \n",
    "\n",
    "        if ID < 56 and par !='cbv':\n",
    "            conditions=['rest', 'control', 'mem', 'calc']\n",
    "\n",
    "        if ID > 55: \n",
    "            conditions=['control', 'calc']\n",
    "\n",
    "        if par == 'cbv' and ID < 56:\n",
    "            conditions=['control']\n",
    "\n",
    "        \n",
    "        for cond in conditions: \n",
    "            if par != 'cmro2':\n",
    "                par_MNI = os.path.join(derivatives_dir, sub, 'qmri',sub + '_task-'+ cond +'_space-MNI152_' + par + '.nii.gz')\n",
    "            if par == 'cmro2' and cond != 'calc':\n",
    "                par_MNI = os.path.join(derivatives_dir, sub, 'qmri',sub + '_task-'+ cond +'_space-MNI152_desc-orig_' + par + '.nii.gz')\n",
    "            if par == 'cmro2' and cond == 'calc':\n",
    "                par_MNI = os.path.join(derivatives_dir, sub, 'qmri',sub + '_task-'+ cond +'_space-MNI152_desc-orig_' + par + '.nii.gz')\n",
    "\n",
    "            T1w_MNIspace = os.path.join(derivatives_dir, sub, 'anat',sub + '_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz')\n",
    "            plotting.plot_img(par_MNI,  bg_img =T1w_MNIspace, cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), threshold = 0.5, black_bg=True, title = par + ' : ' + sub)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### normalize T2map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    ". /etc/fsl/5.0/fsl.sh\n",
    ". /etc/ants/ants.sh\n",
    "\n",
    "##only T2 map\n",
    "parameters=\"T2map\"\n",
    "\n",
    "sids=\"19 20 21 23 26 27 28 29 30 31 32 34 33 35 36 37 38 39 40 43 44 46 47 48 49 50 51 52 54 55 58 59 60 61 63 64 65 66 67 68\"\n",
    "#sids=\"58 59 60 61 63 64 65 66 67 68 \"\n",
    "#sids=\"20\"\n",
    "\n",
    "## T2map, no condition\n",
    "for ID in $sids; \n",
    "do\n",
    "    ##define directories\n",
    "    SID=p0${ID}\n",
    "    echo $SID\n",
    "    base_path=\"/home/tumnic/sepp/neuroenergeticslab_notebooks/erc-wp1/manuscript\" \n",
    "    locfolder=\"${base_path}/data/derivatives/sub-${SID}\"\n",
    "    \n",
    "    ##define input\n",
    "    ## output from fmriprep\n",
    "    T1w=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T1w.nii.gz\"\n",
    "    br_mask=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_brain_mask.nii.gz\"\n",
    "    transform_to_standard=\"${locfolder}/anat/sub-${SID}_from-T1w_to-MNI152NLin6Asym_mode-image_xfm.h5\"\n",
    "    #echo $transform_to_standard\n",
    "    T1w_MNIspace=\"${locfolder}/anat/sub-${SID}_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz\"\n",
    "    #echo $T1w_MNIspace\n",
    "    ## output from qmri\n",
    "    T1w_T2space=\"${locfolder}/qmri/sub-${SID}_space-T2_T1w.nii\"\n",
    "    \n",
    "    ##define output\n",
    "    T1w_brain=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T1w_brain.nii.gz\"\n",
    "    T2_to_T1w=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T2_to_T1w.mat\"\n",
    "    \n",
    "    ##brain-extract T1w\n",
    "    echo \"apply brain mask for subject ${SID} (BET)\"\n",
    "    fslmaths $T1w -mas $br_mask $T1w_brain\n",
    "    #bet $T1w $T1w_brain\n",
    "    \n",
    "    ##get T2 to T1w transformation mat\n",
    "    flirt -in $T1w_T2space -ref $T1w_brain -omat $T2_to_T1w -dof 6\n",
    "\n",
    "    \n",
    "    for par in $parameters;\n",
    "    do    \n",
    "\n",
    "        \n",
    "        parameter_map=\"${locfolder}/qmri/sub-${SID}_space-T2_T2map.nii\"\n",
    "\n",
    "\n",
    "        ##define outputs\n",
    "        out_T1w=\"${locfolder}/qmri/sub-${SID}_space-T1w_T2map.nii.gz\"\n",
    "        out_MNI=\"${locfolder}/qmri/sub-${SID}_space-MNI152_T2map.nii.gz\"\n",
    "\n",
    "\n",
    "\n",
    "           ##apply transformation from lowres T2 to highres T1 space for parameter maps\n",
    "        flirt -in $parameter_map -ref $T1w_brain -applyxfm -init $T2_to_T1w -out $out_T1w\n",
    "\n",
    "        ##apply transformation (saved by fmriprep) from T1 to standard space\n",
    "        antsApplyTransforms -d 3 --float 1 --verbose 1 -i $out_T1w -r $T1w_MNIspace -t $transform_to_standard -o $out_MNI\n",
    "\n",
    "    done\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## normalize semi-quant CMRO2 maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    ". /etc/fsl/5.0/fsl.sh\n",
    ". /etc/ants/ants.sh\n",
    "\n",
    "parameters=\"oef cmro2\"\n",
    "\n",
    "sids=\"19 20 21 23 26 27 28 29 30 31 32 34 33 35 36 37 38 39 40 43 44 46 47 48 49 50 51 52 54 55 58 59 60 61 63 64 65 66 67 68\"\n",
    "#sids=\"19 20 21 23 26 27 28 29 30 31 32 34 33 35 36 37 38 39 40 43 44 46 47 48 49 50 51 52 54 55\"\n",
    "\n",
    "#sids=\"19 20\"\n",
    "\n",
    "for par in $parameters; do\n",
    "for ID in $sids; \n",
    "do\n",
    "    ##define directories\n",
    "    SID=p0${ID}\n",
    "    echo $SID\n",
    "    base_path=\"/home/tumnic/sepp/neuroenergeticslab_notebooks/erc-wp1/manuscript\" \n",
    "    locfolder=\"${base_path}/data/derivatives/sub-${SID}\"\n",
    "    \n",
    "    ##define input\n",
    "    ## output from fmriprep\n",
    "    T1w=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T1w.nii.gz\"\n",
    "    br_mask=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_brain_mask.nii.gz\"\n",
    "    transform_to_standard=\"${locfolder}/anat/sub-${SID}_from-T1w_to-MNI152NLin6Asym_mode-image_xfm.h5\"\n",
    "    #echo $transform_to_standard\n",
    "    T1w_MNIspace=\"${locfolder}/anat/sub-${SID}_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz\"\n",
    "    #echo $T1w_MNIspace\n",
    "    ## output from qmri\n",
    "    T1w_T2space=\"${locfolder}/qmri/sub-${SID}_space-T2_T1w.nii\"\n",
    "    \n",
    "    ##define output\n",
    "    T1w_brain=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T1w_brain.nii.gz\"\n",
    "    T2_to_T1w=\"${locfolder}/anat/sub-${SID}_desc-fmriprep_T2_to_T1w.mat\"\n",
    "    \n",
    "    ##brain-extract T1w\n",
    "    #echo \"apply brain mask for subject ${SID} (BET)\"\n",
    "    #fslmaths $T1w -mas $br_mask $T1w_brain\n",
    "    #bet $T1w $T1w_brain\n",
    "    \n",
    "    ##get T2 to T1w transformation mat\n",
    "    #flirt -in $T1w_T2space -ref $T1w_brain -omat $T2_to_T1w -dof 6\n",
    "\n",
    "    \n",
    "    if [ \"$ID\" -lt \"56\" ]; then\n",
    "        if [ par != \"cbv\" ]; then\n",
    "            #tasks=\"calc mem\"\n",
    "            tasks=\"calc\"\n",
    "            #baseline=\"control\"\n",
    "            baseline=\"rest\"\n",
    "        fi\n",
    "    fi\n",
    "    if [ \"$ID\" -gt \"55\" ]; then\n",
    "            tasks=\"calc\"\n",
    "            baseline=\"control\"\n",
    "    fi\n",
    "\n",
    "        \n",
    "    for task in $tasks; do            \n",
    "                \n",
    "        parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-T2_desc-semi-quant_${par}.nii.gz\"\n",
    "        out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-T1w_desc-semi-quant_${par}.nii.gz\"\n",
    "        out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-MNI152_desc-semi-quant_${par}.nii.gz\"\n",
    "\n",
    "        if [ $task == \"calc\" ]; then\n",
    "                parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-T2_desc-semi-quant_${par}.nii.gz\"\n",
    "                out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-T1w_desc-semi-quant_${par}.nii.gz\"\n",
    "                out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-MNI152_desc-semi-quant_${par}.nii.gz\"\n",
    "        fi\n",
    "            \n",
    "            ##apply transformation from lowres T2 to highres T1 space for parameter maps\n",
    "        flirt -in $parameter_map -ref $T1w_brain -applyxfm -init $T2_to_T1w -out $out_T1w\n",
    "\n",
    "            ##apply transformation (saved by fmriprep) from T1 to standard space\n",
    "        antsApplyTransforms -d 3 --float 1 --verbose 1 -i $out_T1w -r $T1w_MNIspace -t $transform_to_standard -o $out_MNI\n",
    "\n",
    "        ## also normalize CBV-corrected cmro2 maps\n",
    "        if [ $task == \"calc\" ]; then\n",
    "            parameter_map=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-T2_desc-semi-quant-corrected_${par}.nii.gz\"\n",
    "            out_T1w=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-T1w_desc-semi-quant-corrected_${par}.nii.gz\"\n",
    "            out_MNI=\"${locfolder}/qmri/sub-${SID}_task-${task}_base-${baseline}_space-MNI152_desc-semi-quant-corrected_${par}.nii.gz\"\n",
    "            ##apply transformation from lowres T2 to highres T1 space for parameter maps\n",
    "            flirt -in $parameter_map -ref $T1w_brain -applyxfm -init $T2_to_T1w -out $out_T1w\n",
    "\n",
    "            ##apply transformation (saved by fmriprep) from T1 to standard space\n",
    "            antsApplyTransforms -d 3 --float 1 --verbose 1 -i $out_T1w -r $T1w_MNIspace -t $transform_to_standard -o $out_MNI\n",
    "        fi\n",
    "            \n",
    "\n",
    "        \n",
    "    done\n",
    "done\n",
    "done\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## show normalized parameter maps\n",
    "parameters = ['cmro2', 'oef']\n",
    "\n",
    "for par in parameters:\n",
    "    for ID in sids:\n",
    "        sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "        print(sub)\n",
    "        sub_ID = \"p{:03d}\".format(ID) \n",
    "\n",
    "        if ID < 56 and par !='cbv':\n",
    "            tasks=['mem', 'calc']\n",
    "            baseline='control'\n",
    "\n",
    "        if ID > 55: \n",
    "            tasks=[ 'calc']\n",
    "            baseline='control'\n",
    "\n",
    "\n",
    "\n",
    "        for cond in tasks:\n",
    "            par_MNI = os.path.join(derivatives_dir, sub, 'qmri',sub + '_task-'+ cond +'_base-'+baseline+'_space-MNI152_desc-semi-quant_' + par + '.nii.gz')\n",
    "            \n",
    "            T1w_MNIspace = os.path.join(derivatives_dir, sub, 'anat',sub + '_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz')\n",
    "            plotting.plot_img(par_MNI,  bg_img =T1w_MNIspace, cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), threshold = 0, black_bg=True, title = par + ' ' + sub)\n",
    "        \n",
    "        \n",
    "            if cond == 'calc':\n",
    "                par_MNI = os.path.join(derivatives_dir, sub, 'qmri',sub + '_task-'+ cond +'_base-'+baseline+'_space-MNI152_desc-semi-quant-corrected_' + par + '.nii.gz')\n",
    "\n",
    "                T1w_MNIspace = os.path.join(derivatives_dir, sub, 'anat',sub + '_space-MNI152NLin6Asym_res-2_desc-preproc_T1w.nii.gz')\n",
    "                plotting.plot_img(par_MNI,  bg_img =T1w_MNIspace, cmap='jet',display_mode='z', cut_coords = (0, 10, 20, 30, 40, 50), threshold = 0, black_bg=True, title = par + ' corrected ' + sub)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## threshold native-space parameter maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## run qBOLD toolbox\n",
    "%run qBOLD_fun.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#grey matter mask\n",
    "GMthresh = 0.5; #mandatory\n",
    "\n",
    "#mask CSF influenced areas\n",
    "T2thresh = 90; #mandatory\n",
    "\n",
    "#mask R2strich artifact areas\n",
    "R2thresh = 9 ; #mandatory\n",
    "\n",
    "#mask CBV \n",
    "CBVthresh = 10; #optional: comment out if not used, default = 20\n",
    "\n",
    "#mask T2S\n",
    "T2Sthresh = 90; #optional: comment out if not used, default = 150\n",
    "\n",
    "#mask OEF\n",
    "OEFthresh = 0.9; #optional: comment out if not used, default = 1\n",
    "\n",
    "#mask CBF\n",
    "CBFthresh = 90; #optional: comment out if not used, default = 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters=['T2map', 'T2Smap', 'cbv', 'cbf', 'R2prime', 'oef', 'cmro2']\n",
    "\n",
    "for p, par in enumerate(parameters):\n",
    "    for i, ID in enumerate(sids): #loop over subjects\n",
    "        sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "        print(sub)\n",
    "\n",
    "        dir_deriv = os.path.join(data_dir, 'derivatives', sub)\n",
    "    \n",
    "\n",
    "\n",
    "        if ID < 56:\n",
    "            conds = ['rest', 'control', 'mem', 'calc']\n",
    "        if ID > 55:\n",
    "            conds = ['control', 'calc']   \n",
    "\n",
    "        mask_list = create_qBOLD_masks(sub, dir_deriv, conds, GMthresh, T2thresh, R2thresh, CBVthresh, T2Sthresh, OEFthresh,  CBFthresh)\n",
    "        \n",
    "        if par == 'T2map':\n",
    "            conds = ['control']\n",
    "        if par == 'func':\n",
    "            conds = ['all']\n",
    "        if par == 'cbv' and ID < 56:\n",
    "            conds=['control']\n",
    "        if par == 'cbv' and ID > 55:\n",
    "            conds=['control', 'calc']\n",
    "            \n",
    "\n",
    "        for cond in conds:\n",
    "            ## load parameter map\n",
    "            if par != 'T2map' and par != 'cmro2':\n",
    "                    par_nii =  os.path.join(dir_deriv, 'qmri', sub + '_task-'+cond+'_space-T2_'+ par +'.nii.gz')\n",
    "            if par == 'T2map':\n",
    "                    par_nii =  os.path.join(dir_deriv, 'qmri', sub + '_space-T2_'+ par +'.nii.gz')\n",
    "            if par == 'cmro2' and cond != 'calc':\n",
    "                    par_nii =  os.path.join(dir_deriv, 'qmri', sub + '_task-'+cond+'_space-T2_desc-orig_'+ par +'.nii.gz')\n",
    "            if par == 'cmro2' and cond == 'calc':\n",
    "                    par_nii =  os.path.join(dir_deriv, 'qmri', sub + '_task-'+cond+'_space-T2_desc-CBV_'+ par +'.nii.gz')\n",
    "            if par == 'func':\n",
    "                    par_nii =  os.path.join(dir_deriv, 'func', sub + '_task-'+cond+'_space-T2_filtered_func.nii.gz')\n",
    "                    \n",
    "            ##mask parameter map\n",
    "            out_name = par_nii[:-7] + '_qBmasked.nii.gz'\n",
    "            par_masked = mask_nii(par_nii, mask_list, out_name)\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## mask 1st level ROIs, output from 1st level GLM\n",
    "\n",
    "for i, ID in enumerate(sids): #loop over subjects\n",
    "    sub = \"sub-p{:03d}\".format(ID) # subject id, eg 'p021'\n",
    "    print(sub)\n",
    "\n",
    "    dir_deriv = os.path.join(derivatives_dir, sub)\n",
    "            \n",
    "    if ID<56:\n",
    "        ROIs = ['calccontrol', 'controlcalc', 'calcrest','restcalc', 'memcontrol', 'controlmem'] #visuomotor: control>rest\n",
    "    if ID>55:\n",
    "        ROIs = ['calccontrol', 'controlcalc']\n",
    "\n",
    "    mask_list = [ os.path.join(dir_deriv, 'qmri',  sub + '_space-T2_GM_mask.nii.gz'), \n",
    "             os.path.join(dir_deriv, 'qmri',  sub + '_space-T2_T2map_mask.nii.gz'), \n",
    "             os.path.join(dir_deriv, 'qmri',  sub + '_space-T2_R2prime_mask.nii.gz'),\n",
    "            os.path.join(dir_deriv, 'qmri', sub + '_task-control_space-T2_cbv.nii.gz'),\n",
    "             os.path.join(dir_deriv, 'qmri',  sub + '_space-T2_T2Smap_mask.nii.gz'), \n",
    "             os.path.join(dir_deriv, 'qmri',  sub + '_space-T2_oef_mask.nii.gz'), \n",
    "              os.path.join(dir_deriv, 'qmri',  sub + '_space-T2_cbf_mask.nii.gz')]\n",
    "\n",
    "\n",
    "    for ROI in ROIs:\n",
    "        ROI_nii = os.path.join(dir_deriv, 'func', sub + '_1stlevel_'+ROI+'_space-T2.nii.gz')\n",
    "                    \n",
    "        ##mask parameter map\n",
    "        out_name = ROI_nii[:-7] + '_qBmasked.nii.gz'\n",
    "        par_masked = mask_nii(ROI_nii, mask_list, out_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PR postPLS",
   "language": "python",
   "name": "pr_postpls_hist"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}