notebooks/extract_all_timeseries.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "99aae7af",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from nilearn import datasets, image, plotting, masking, maskers\n",
    "from bids import BIDSLayout\n",
    "import json\n",
    "from io import StringIO\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "\n",
    "# from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms\n",
    "# from niworkflows.interfaces.itk import MultiApplyTransforms\n",
    "# from nilearn.connectome import ConnectivityMeasure\n",
    "# from sklearn.base import clone, BaseEstimator, TransformerMixin\n",
    "# from sklearn.model_selection import StratifiedShuffleSplit, StratifiedGroupKFold\n",
    "# from sklearn.linear_model import LogisticRegression\n",
    "# from sklearn.feature_selection import SelectKBest\n",
    "# from sklearn.decomposition import TruncatedSVD\n",
    "# from sklearn.compose import ColumnTransformer\n",
    "# from sklearn.pipeline import Pipeline\n",
    "# from sklearn.preprocessing import FunctionTransformer\n",
    "# import sklearn.datasets\n",
    "# from sklearn import svm\n",
    "# from sklearn.metrics import accuracy_score\n",
    "# from scipy.special import btdtr\n",
    "# from statsmodels.stats.multitest import multipletests\n",
    "# import seaborn as sns\n",
    "# from matplotlib import pyplot as plt\n",
    "# import itertools\n",
    "# import dask_ml.model_selection\n",
    "# import sklearn.model_selection\n",
    "\n",
    "from joblib import Memory\n",
    "import os\n",
    "from joblib import Parallel, delayed\n",
    "%matplotlib inline\n",
    "pd.set_option('display.max_rows', 500)\n",
    "pd.set_option('display.max_columns', 500)\n",
    "pd.set_option('display.max_colwidth', 500)\n",
    "pd.set_option('display.width', 1000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9afd800b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_list_merge(df, ll, col_name):\n",
    "    dfs = []\n",
    "    for lx in ll:\n",
    "        tmp = df.copy()\n",
    "        tmp[col_name] = lx\n",
    "        dfs.append(tmp)\n",
    "    return pd.concat(dfs).reset_index(drop=True)\n",
    "\n",
    "def select_confounds(cfds_path, cfds_sel):\n",
    "    cfd = pd.read_csv(cfds_path, sep='\\t')\n",
    "    cols = []\n",
    "    #cfd = cfd.drop('censored', axis=1)\n",
    "    if \"-censor\" in cfds_sel: cols.extend([cc for cc in cfd.columns if 'censor' in cc])\n",
    "    if \"-cosine\" in cfds_sel: cols.extend([cc for cc in cfd.columns if 'cosine' in cc])\n",
    "    if \"-aroma\" in cfds_sel: cols.extend([cc for cc in cfd.columns if 'aroma' in cc])\n",
    "    if \"-motion\" in cfds_sel: cols.extend([cc for cc in cfd.columns if ('trans' in cc) or ('rot' in cc)])\n",
    "    if \"-physio\" in cfds_sel: cols.extend([cc for cc in cfd.columns if 'phys' in cc])\n",
    "    cfds = cfd.loc[:, cols].copy()\n",
    "    return cfds.values\n",
    "\n",
    "def make_path(ents, updates, pattern, derivatives_dir, build_path, check_exist=True, derivatives=False,\n",
    "              check_parent=True, mkdir=False, make_parent=False):\n",
    "    mp_ents = ents.copy()\n",
    "    if derivatives:\n",
    "        mp_dir = fmriprep_out / f'sub-{ents[\"subject\"]}/out'\n",
    "    else:\n",
    "        mp_dir = derivatives_dir \n",
    "    mp_ents.update(updates)\n",
    "    mp_file = mp_dir / build_path(mp_ents, pattern, validate=False, absolute_paths=False)\n",
    "    if check_exist and not mp_file.exists():\n",
    "        raise FileNotFoundError(mp_file.as_posix())\n",
    "    elif check_parent and not mp_file.parent.exists():\n",
    "         raise FileNotFoundError(mp_file.parent.as_posix())\n",
    "    if mkdir:\n",
    "        mp_file.mkdir(parents=True, exist_ok=True)\n",
    "    elif make_parent:\n",
    "        mp_file.parent.mkdir(parents=True, exist_ok=True)\n",
    "    return mp_file\n",
    "\n",
    "def r_to_p(r, n):\n",
    "    from scipy.special import btdtr\n",
    "    \"\"\"\n",
    "    Adapted from Scipy Stats: https://github.com/scipy/scipy/blob/38cd478334a07b6e9dfe224a3ff66ad356c53524/scipy/stats/stats.py#L4069-L4070\n",
    "    \"\"\"\n",
    "    r = np.array(r).copy()\n",
    "    ab = n / 2 - 1\n",
    "    prob = 2*btdtr(ab, ab, 0.5 * (1 - np.abs(np.float64(r))))\n",
    "    return prob\n",
    "\n",
    "def test_r_to_p():\n",
    "    rng = np.random.default_rng()\n",
    "    n = 100\n",
    "    x = rng.normal(size=n)\n",
    "    y = rng.normal(size=n)\n",
    "    r, true_p = stats.pearsonr(x,y)\n",
    "    assert true_p == r_to_p(r, n)\n",
    "    y = rng.normal(size=n) + x * 0.3\n",
    "    r, true_p = stats.pearsonr(x,y)\n",
    "    assert true_p == r_to_p(r, n)\n",
    "test_r_to_p()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "db83fc8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "project_root = Path('..')\n",
    "bids_dir = project_root / 'data/CATD'\n",
    "derivatives_dir = bids_dir / 'derivatives'\n",
    "fmriprep_out = derivatives_dir / 'fmriprep' / 'fmriprep_v21.0.0'\n",
    "\n",
    "nilearn_out = derivatives_dir / 'rest_processed'\n",
    "nilearn_out.mkdir(exist_ok=True)\n",
    "jobids = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d979c0fd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.21 ms, sys: 2.18 ms, total: 4.4 ms\n",
      "Wall time: 4.54 ms\n"
     ]
    }
   ],
   "source": [
    "prep_df = pd.read_csv(derivatives_dir / 'summary_tables' / 'prep_df.csv')\n",
    "rest_df = prep_df.loc[prep_df.task == 'rest'].copy()\n",
    "rest_paths = pd.read_csv(derivatives_dir / 'summary_tables' / 'rest_sessamplenothresh_prep.csv')\n",
    "rest_paths['entities'] = rest_paths['entities'].apply(lambda x: json.loads(x.replace(\"'\", '\"')))\n",
    "\n",
    "database_path='./pybids140_public'\n",
    "%time layout = BIDSLayout(bids_dir, database_path=database_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3c2d0849-e004-4a86-a153-adb8bd75ed36",
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_paths['nv_noncensored'] = (rest_paths.nv - 4) * (1-rest_paths.pct_censored)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5946806b-b23e-4743-81ee-62be39ef49eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_paths['above_mfq_thresh'] = rest_paths.s_mfq_tot >= 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ce89dd89-67cf-4613-a8c2-dded7d62765a",
   "metadata": {},
   "outputs": [],
   "source": [
    "path_cols = ['file'] + list(rest_paths.columns[rest_paths.columns.str.contains('path')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3da3a19e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# rearrange paths to be longish\n",
    "keep_columns_notd = [\n",
    "    'subject',\n",
    "    'run',\n",
    "    'suffix',\n",
    "    'session',\n",
    "    'task',\n",
    "    'space',\n",
    "    'datatype',\n",
    "    'extension',\n",
    "    'has_phys',\n",
    "    'sex',\n",
    "    'group',\n",
    "    'age',\n",
    "    'scanner',\n",
    "    'inpatient',\n",
    "    'acq_date',\n",
    "    'usable',\n",
    "    'dropout',\n",
    "    'nv_noncensored',\n",
    "    'above_mfq_thresh',\n",
    "    'entities',\n",
    "    'file',\n",
    "    'mask_path',\n",
    "    'cfd_new_path',\n",
    "#     'rest_corr_path',\n",
    "#     'rest_ts_path'\n",
    "]\n",
    "rest_notd = rest_paths.loc[:, keep_columns_notd]\n",
    "rest_notd['processing'] = 'notd'\n",
    "# rest_notd = rest_notd.rename(columns={'rest_corr_path': 'orig_rest_corr_path'})\n",
    "# rest_notd = rest_notd.rename(columns={'rest_ts_path': 'orig_rest_ts_path'})\n",
    "\n",
    "rest_p = rest_notd.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bb9d60dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# atlasses = [image.load_img((nilearn_out / 'basc122_with_schaefer400_ldlpfc.nii.gz').as_posix()),\n",
    "#             image.load_img((nilearn_out / 'Schaefer2018_400Parcels_7Networks_order_FSLMNI152_2mm.nii.gz').as_posix())]\n",
    "# atlasses = {'Basc122': image.load_img((nilearn_out / 'basc122_with_schaefer400_ldlpfc.nii.gz').as_posix()),\n",
    "#             'Schaefer400': image.load_img((nilearn_out / 'Schaefer2018_400Parcels_7Networks_order_FSLMNI152_2mm.nii.gz').as_posix())}\n",
    "atlasses = {'Basc122': (nilearn_out / 'basc122_flat_with_ldlpfc.nii.gz').as_posix(),\n",
    "            }\n",
    "rest_pa = df_list_merge(rest_p, atlasses.keys(), 'atlas')\n",
    "rest_pa['atlas_img'] = rest_pa.atlas.replace(atlasses)\n",
    "confounds = [\"-censor\", \"-cosine\", \"-motion\", \"-physio\"]\n",
    "rest_pa['counfound_array'] = rest_pa.cfd_new_path.apply(lambda x: select_confounds(x, confounds))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7533c0d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_pa['desc'] = rest_pa.processing+rest_pa.atlas\n",
    "rest_pa['pa_ents'] = rest_pa.loc[:, ['subject', 'session', 'task', 'run', 'space', 'suffix', 'datatype', 'desc']].to_dict('records')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "feac8496",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_run_name = 'ts_extract_allmodels'\n",
    "cfds_sel=['censor','cosine','motion','physio']\n",
    "ts_out = nilearn_out / ts_run_name\n",
    "ts_out.mkdir(exist_ok=True)\n",
    "ts_pattern = 'sub-{subject}/ses-{session}/sub-{subject}_ses-{session}_task-{task}_run-{run}_space-{space}_desc-{desc}_{suffix}.{extension}'\n",
    "ts_updates = {'suffix':'timeseries', 'extension':'csv'}\n",
    "rest_pa['rest_ts_path'] = rest_pa.pa_ents.apply(lambda x: make_path(x,\n",
    "                                                            ts_updates,\n",
    "                                                            ts_pattern,\n",
    "                                                            ts_out,\n",
    "                                                            layout.build_path,\n",
    "                                                            check_exist=False,\n",
    "                                                            check_parent=False,\n",
    "                                                            make_parent=True\n",
    "                                                                          ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "40111b7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_ts_and_save(row, n_dummy=4):\n",
    "    regions_extracted_img = image.load_img(row.atlas_img)\n",
    "    img = image.load_img(row.file)\n",
    "    \n",
    "    masker = maskers.NiftiLabelsMasker(regions_extracted_img,\n",
    "                                          mask_img=row.mask_path,\n",
    "                                          resampling_target=\"data\",\n",
    "                                          low_pass=.1,\n",
    "                                          high_pass=.01,\n",
    "                                          t_r=2.5)\n",
    "    extracted = masker.fit_transform(row.file, confounds=row.counfound_array)[n_dummy:,:]\n",
    "    pd.DataFrame(extracted).to_csv(row.rest_ts_path, header=None, index=None)\n",
    "    return extracted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "18bfcccb-cb5e-42e6-9263-f4e00cfcfa9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_exist = rest_pa.rest_ts_path.apply(lambda x: x.exists())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1d4d7dd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "extract_jobs = rest_pa.loc[~ts_exist].apply(lambda row: delayed(extract_ts_and_save)(row.copy()), axis=1)\n",
    "#extract_jobs = rest_pa.apply(lambda row: delayed(extract_ts_and_save)(row.copy()), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1c82ec66-656e-42ab-9da7-0373d98171e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "504"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(extract_jobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9a4fd188-c71e-4f44-817e-35dfacd8381e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n",
      "/gpfs/gsfs11/users/MBDU/mid_share/env/lib/python3.9/site-packages/scipy/ndimage/measurements.py:684: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return sum / numpy.asanyarray(count).astype(numpy.float64)\n"
     ]
    }
   ],
   "source": [
    "if len(extract_jobs) > 0:\n",
    "    res = Parallel(n_jobs=20)(extract_jobs)\n",
    "extracted = rest_pa.loc[ts_exist, 'rest_ts_path'].apply(lambda x: pd.read_csv(x,header=None).values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "88c38aee",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# make sure extraction worked\n",
    "ts_exist = rest_pa.rest_ts_path.apply(lambda x: x.exists())\n",
    "assert ts_exist.all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3c0d7995",
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_pa['rest_ts'] = extracted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "fff1a547-eb2f-4aa1-97d0-fe2819a29351",
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_pa.to_csv(derivatives_dir / 'summary_tables' / 'rest_paths_with_timeseries.csv', index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b452157",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}