1_make_dataframe.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "96007540",
   "metadata": {},
   "source": [
    "# Making the dataframe\n",
    "## This must be run before other code sections\n",
    "## This code contains\n",
    "1) Phenotypic data loading\n",
    "2) Quality Control\n",
    "3) Calculating covariates and FA metrics\n",
    "\n",
    "### If replicating this study, make sure your dataset is BIDS-compliant and:\n",
    "1) Update the variable \"bids_dir\" to reflect where your data is stored\n",
    "2) Have your phenotypic data stored in the BIDS code directory with the name \"HBN_query.csv\"\n",
    "3) Make sure you have collected at the very least: Identifiers, Basic demographics, Clinican Diagnoses, TOWRE scores, and EHQ scores from the HBN data portal\n",
    "\n",
    "### This could be adapted for other studies that use this suite of neuroimaging scripts on different subjects:\n",
    "1) Make sure your phenotypic file is organized in the same way (one row per subject, one column per trait)\n",
    "2) Either change the headings of your phenotypic file to match the HBN headers, or change the script to match your headers\n",
    "3) Point the variable \"HBN_query\" to your phenotypic file\n",
    "4) Change the \"representative_subject\" variable to one with good quality data and a complete gradient table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "811b9ae4",
   "metadata": {},
   "source": [
    "# 1) Import packages and load data\n",
    "### Make sure to read the instructions in the github and the preamble to this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "566484ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import glob\n",
    "import filecmp\n",
    "import json\n",
    "import nilearn\n",
    "import nilearn.image\n",
    "import fsl\n",
    "import fsl.wrappers\n",
    "import fsl.utils.image.resample\n",
    "from fsl.data.image import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44018e21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CHANGE THIS VARIABLE TO REFLECT WHERE YOUR HBN BIDS DATA LIVE ()\n",
    "bids_dir = '/' # Path should end with a '/'\n",
    "\n",
    "# IF USING YOUR OWN DATASET, CHANGE THIS VARIABLE TO A SUBJECT WITH A GOOD QUALITY RUN AND COMPLETE GRADIENT TABLE\n",
    "representative_subject = 'sub-NDARAA948VFH/'\n",
    "num_bad_slices_thresh = 9 # Any subjects with more bad slices than this are excluded (see QC). Change this if you wish\n",
    "\n",
    "# TRACTS TO ANALYZE\n",
    "tracts = ['AF_left','AF_right','SLF_I_left','SLF_I_right','SLF_II_left','SLF_II_right',\n",
    "         'SLF_III_left','SLF_III_right','ILF_left','ILF_right','IFO_left','IFO_right',\n",
    "         'UF_left','UF_right','SCP_left','SCP_right','ICP_left','ICP_right','MCP','CC_7']\n",
    "\n",
    "# THESE VARIABLES SHOULD NOT CHANGE IF THE ANALYSIS WAS RUN ACCORDING TO INSTRUCTIONS\n",
    "code_dir = bids_dir+'code/'\n",
    "derivatives_dir = bids_dir+'derivatives/'\n",
    "qsiprep_dir = derivatives_dir+'qsiprep/'\n",
    "tractseg_dir = derivatives_dir+'TractSeg_NEW/'\n",
    "freesurfer_dir = derivatives_dir+'freesurfer/'\n",
    "\n",
    "# LOAD THE PHENOTYPIC DATA\n",
    "HBN_query = pd.read_csv(code_dir+'HBN_query2.csv')\n",
    "\n",
    "# GET SUBJECTS WITH PHENOTYPIC DATA\n",
    "query_subs = ['sub-'+name[:-11] for name in HBN_query['Identifiers']]\n",
    "\n",
    "# GET SUBJECTS WITH TractSeg SEGMENTATION DATA\n",
    "tractseg_subs = np.transpose([path.split('/')[-1] for path in glob.glob(tractseg_dir+'sub-*')])\n",
    "tractseg_subs.sort()\n",
    "\n",
    "# GET SUBJECTS WITH FreeSurfer STATS DATA\n",
    "freesurfer_subs = [string.split('/')[-3] for string in glob.glob(freesurfer_dir+'sub-*/stats/aseg.stats')]\n",
    "freesurfer_subs.sort()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e5a6342",
   "metadata": {},
   "source": [
    "# 2) Find intersection of subjects with valid TOWRE scores, handedness data, FreeSurfer, and segmentation data\n",
    "\n",
    "### This assumes any folder in your TractSeg directory beginning with \"sub-\" has complete data. Please remove failed TractSeg subjects from your TractSeg directory or rename them (e.g. \"BAD_sub-XX\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7e7354d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GET SUBJECTS WITH VALID TOWRE SCORES, HANDEDNESS DATA\n",
    "towre_query_subs = np.asarray(query_subs)[HBN_query['TOWRE,TOWRE_Valid']=='1']\n",
    "hand_query_subs = np.asarray(query_subs)[HBN_query['EHQ,Data_entry']=='Complete']\n",
    "\n",
    "# GET INTERSECTION OF SUBJECTS WITH PHENOTYPIC AND IMAGING DATA\n",
    "subs_pheno_img = list(set(tractseg_subs) & set(freesurfer_subs) & set(towre_query_subs) & set(hand_query_subs))\n",
    "subs_pheno_img.sort()\n",
    "\n",
    "# FIND WHERE THESE SUBJECTS ARE LOCATED IN THE PHENOTYPIC FILE\n",
    "pheno_img_query_indexes = np.transpose([np.where(np.asarray(query_subs)==sub)[0][0] for sub in subs_pheno_img])\n",
    "\n",
    "# EXTRACT THE TOWRE SCORES\n",
    "towre_scores_with_nan=HBN_query['TOWRE,TOWRE_Total_Scaled'][pheno_img_query_indexes]\n",
    "\n",
    "# FILTER OUT THE NANs\n",
    "good_inds=np.asarray(towre_scores_with_nan.isnull()==False)\n",
    "towre_scores = np.asarray(list(map(int, np.asarray(towre_scores_with_nan[good_inds]))))\n",
    "\n",
    "# GET LIST OF SUBEJCTS WITH AND FULL PHENOTYPIC DATA AND IMAGING DATA\n",
    "subs_pheno_img = np.asarray(subs_pheno_img)[good_inds]\n",
    "\n",
    "# GET THEIR LOCATIONS IN THE PHENOTYPIC FILE\n",
    "pheno_img_query_indexes = pheno_img_query_indexes[good_inds]\n",
    "print('A total of',np.size(subs_pheno_img),'subjects have full phenotypic and segmentation data.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c78a002",
   "metadata": {},
   "source": [
    "# 3) Quality Control\n",
    "### Here we use the QC outputs of QSIPrep to exclude any subjects with poor quality images, inconsistent diffusion parameters, or an age > 18.\n",
    "### If you have additional subjects you would like to exclude (e.g. visually analyzing QSIPrep HTML files or TractSeg outputs), add their full subject name (e.g. sub-PROJECTXXXX) to the variable \"more_subs_exclude\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1848a756",
   "metadata": {},
   "outputs": [],
   "source": [
    "more_subs_exclude = []\n",
    "\n",
    "# INITIALIZE OUTPUT VARIABLES\n",
    "subs_qc = [] # sub names\n",
    "inds_qc = [] # phenotypic file indexes\n",
    "\n",
    "# LOAD QSIPREP GROUP QC OUTPUT\n",
    "qc_file  = qsiprep_dir+'dwiqc.json'\n",
    "with open(qc_file) as f:\n",
    "        qc_data = json.load(f)['subjects']\n",
    "# GET ORDER OF SUBJECTS in QC Data\n",
    "qc_sub_order = np.asarray([sub_data['subject_id'] for sub_data in qc_data])\n",
    "# LOAD REPRESENTATIVE B-VALUES TO COMPARE AGAINST\n",
    "b_val = glob.glob(qsiprep_dir+representative_subject+'/dwi/*.bval')[0]\n",
    "\n",
    "# RUN QC\n",
    "print('SUBJECTS BELOW HAVE BEEN EXCLUDED, THE THREE CRITERIA ARE LISTED NEXT TO THE SUBJECT NAME')\n",
    "for ind in range(np.size(subs_pheno_img)):\n",
    "    \n",
    "    # EXTRACT SUBJECT NAME AND INDEX\n",
    "    sub = subs_pheno_img[ind]\n",
    "    query_ind = pheno_img_query_indexes[ind]\n",
    "    \n",
    "    # EXCLUDE SUBJECT IF MANUALLY ADDED TO EXCLUSION LIST\n",
    "    if sub in more_subs_exclude:\n",
    "        continue\n",
    "        \n",
    "    # COMPARE B-VALUES TO REPRESENTATIVE SUBJECT\n",
    "    b_val_sub = glob.glob(qsiprep_dir+sub+'/dwi/*.bval')[0]\n",
    "    b_val_comparison = filecmp.cmp(b_val,b_val_sub)\n",
    "    \n",
    "    # LOAD QC DATA AND CHECK NUMBER OF BAD SLICES\n",
    "    qc_ind = [index for index in range(np.size(qc_sub_order)) if qc_sub_order[index] == sub][0]\n",
    "    num_bad_slices = qc_data[qc_ind]['t1_num_bad_slices']\n",
    "    \n",
    "    # CHECK AGE\n",
    "    age_check = np.asarray(HBN_query['Basic_Demos,Age'][query_ind]) < 19\n",
    "    \n",
    "    # IF QC PASSSES, SAVE OUT SUBJECT\n",
    "    if num_bad_slices <= num_bad_slices_thresh and b_val_comparison and age_check:\n",
    "        subs_qc.append(sub)\n",
    "        inds_qc.append(pheno_img_query_indexes[ind])\n",
    "    # IF NOT, PRINT SUBJECT NAME AND WHY THEY DID NOT PASS\n",
    "    else: \n",
    "        print(sub, 'num_bad_slices =',num_bad_slices, ', b_val_comparison =',b_val_comparison, ', age_check =',age_check)\n",
    "        \n",
    "# CONVERT TO A NUMPY ARRAY FOR CONVENIENCE\n",
    "subs_qc = np.asarray(subs_qc)\n",
    "inds_qc = np.asarray(inds_qc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb30aa3c",
   "metadata": {},
   "source": [
    "# 4) Compute extra covariates\n",
    "### These include total brain volume from FreeSurfer, as well as global FA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "281d3847",
   "metadata": {},
   "outputs": [],
   "source": [
    "# UNCOMMENT TO LOAD PREVIOUSLY CALCULATED VALUES\n",
    "icv = np.load('icv.npy')\n",
    "tbv = np.load('tbv.npy')\n",
    "gFA = np.load('gFA.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93ea66c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "icv = [] # Intracranial volume\n",
    "tbv = [] # Total brain volume\n",
    "gFA = [] # Global FA\n",
    "for sub in subs_qc:\n",
    "    # Load FreeSurfer Stats\n",
    "    fs_stats_path = freesurfer_dir+sub+'/stats/aseg.stats'\n",
    "    with open(fs_stats_path) as f:\n",
    "        lines = f.readlines()\n",
    "    # Extract ICV and add to list\n",
    "    icv_text = lines[34]\n",
    "    icv_sub = float((icv_text.split(',')[-2]))\n",
    "    icv.append(icv_sub)\n",
    "     # Extract total brain volume and add to list\n",
    "    tbv_text = lines[13]\n",
    "    tbv_sub = float((tbv_text.split(',')[-2]))\n",
    "    tbv.append(tbv_sub)\n",
    "\n",
    "    # Calculate global FA\n",
    "    # define image paths  \n",
    "    fa_image_path = tractseg_dir+sub+'/dwi/'+sub+'_space-reorient_desc-preproc_FA.nii.gz'\n",
    "    wm_prob_path = qsiprep_dir+sub+'/anat/'+sub+'_label-WM_probseg.nii.gz'\n",
    "    wm_prob_path = qsiprep_dir+sub+'/anat/'+sub+'_label-WM_probseg.nii.gz'\n",
    "    wm_prob_path_reorient = tractseg_dir+sub+'/dwi/'+sub+'_space-reorient_label-WM_probseg.nii.gz'\n",
    "    \n",
    "    # call fslreorient2std command line from notebook to reorient white matter mask\n",
    "    cmd = 'fslreorient2std '+wm_prob_path+ ' '+ wm_prob_path_reorient\n",
    "    os.system(cmd)\n",
    "    \n",
    "    # this white matter mask file will be created\n",
    "    wm_mask_path = wm_prob_path_reorient.replace('probseg','mask')\n",
    "\n",
    "    # Binarize the probseg map to get a white matter mask, and load that image\n",
    "    fsl.wrappers.fslmaths(wm_prob_path_reorient).thr(0.5).bin().run(wm_mask_path)\n",
    "    wm_mask_img = Image(wm_mask_path)\n",
    "    \n",
    "    # Load the FA, resample it to the mask\n",
    "    fa_img = Image(fa_image_path)\n",
    "    fa_resampled = fsl.utils.image.resample.resampleToReference(fa_img, wm_mask_img)[0]\n",
    "    \n",
    "    # Find average FA within the white matter mask\n",
    "    wm_mask_inds = wm_mask_img.data==1\n",
    "    fa_wm_vals = fa_resampled[wm_mask_inds]\n",
    "    fa_global = np.nanmean(fa_wm_vals)\n",
    "    gFA.append(fa_global)\n",
    "    \n",
    "# Convert to numpy arrays and save out    \n",
    "icv = np.array(icv)\n",
    "tbv = np.array(tbv)\n",
    "gFA = np.array(gFA)\n",
    "np.save('icv.npy', icv)\n",
    "np.save('tbv.npy', tbv)\n",
    "np.save('gFA.npy', gFA)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd3f10e5",
   "metadata": {},
   "source": [
    "# 5) Load demographic data only for QC subjects\n",
    "### To distinguish between typical and poor readers, we use the conventional diagnostic criteria of an age-standardized TOWRE score threshold of 85"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7f3e554",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GET PHENOTYPE DATA OF QCed SUBJECTS\n",
    "subs = subs_qc\n",
    "sex = np.asarray(HBN_query['Basic_Demos,Sex'][inds_qc])\n",
    "age = np.asarray(HBN_query['Basic_Demos,Age'][inds_qc])\n",
    "hand = np.asarray(HBN_query['EHQ,EHQ_Total'][inds_qc]).astype(float)\n",
    "towre = np.asarray(HBN_query['TOWRE,TOWRE_Total_Scaled'][inds_qc]).astype(float)\n",
    "swe = np.asarray(HBN_query['TOWRE,TOWRE_SWE_Scaled'][inds_qc]).astype(float)\n",
    "pde = np.asarray(HBN_query['TOWRE,TOWRE_PDE_Scaled'][inds_qc]).astype(float)\n",
    "n = np.size(subs_qc)\n",
    "\n",
    "# GET CLINICIAN DIAGNOSES\n",
    "inds_tr = []\n",
    "inds_rd = []\n",
    "for index in range(n):\n",
    "    ind_query = inds_qc[index]\n",
    "    dxs = [HBN_query[key][ind_query] for key in HBN_query.keys() if 'DX' in key]\n",
    "    dx_check = [('Impairment in Reading' in dx) for dx in dxs if type(dx)==str]\n",
    "    if sum(dx_check)>0:\n",
    "        inds_rd.append(True)\n",
    "        inds_tr.append(False)         \n",
    "    else:\n",
    "        inds_tr.append(True)\n",
    "        inds_rd.append(False)\n",
    "\n",
    "inds_tr = np.asarray(inds_tr)\n",
    "inds_rd = np.asarray(inds_rd)\n",
    "\n",
    "# SPLIT INTO HIGH AND LOW T\n",
    "towre_thresh = 85 # Conventional diagnostic criterion\n",
    "inds_high_t = (towre >= towre_thresh)\n",
    "inds_low_t = (towre < towre_thresh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "054a9152",
   "metadata": {},
   "source": [
    "# 6) Compute average FA values in tracts\n",
    "## Method 1 (used in paper) - Average within tract segmentation mask\n",
    "## Method 2 - Average tract profiles derived from streamlines (e.g. across 100 nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "998eda0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# UNCOMMENT TO LOAD PREVIOUSLY CALCULATED VALUES\n",
    "tract_vals = np.load('tract_vals.npy',allow_pickle=True)[()]\n",
    "tract_vals2 = np.load('tract_vals2.npy',allow_pickle=True)[()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5057deac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# INITIALIZE OUTPUT DICTIONARIES\n",
    "tract_vals={}\n",
    "tract_vals2={}\n",
    "for tract in tracts:\n",
    "    tract_vals[tract] = []\n",
    "    tract_vals2[tract] = []\n",
    "\n",
    "for sub in subs:\n",
    "    # Load the FA (METHOD 1)\n",
    "    fa_image_path = tractseg_dir+sub+'/dwi/'+sub+'_space-reorient_desc-preproc_FA.nii.gz'\n",
    "    fa_vals = np.asarray(nilearn.image.load_img(fa_image_path).dataobj)\n",
    "    \n",
    "    # Load the FA tract profiles (METHOD 2)\n",
    "    fa_tractometry_path = tractseg_dir+sub+'/FA_tractometry.csv'\n",
    "    fa_tractometry_vals = pd.read_csv(fa_tractometry_path,delimiter=';')\n",
    "    \n",
    "    for tract in tracts:\n",
    "        # METHOD 1\n",
    "        # Load the tract mask\n",
    "        tract_mask_path = tractseg_dir+sub+'/bundle_segmentations/'+tract+'.nii.gz'\n",
    "        tract_mask_vals = np.asarray(nilearn.image.load_img(tract_mask_path).dataobj)\n",
    "        # Isolate FA values within the mask\n",
    "        tract_mask_inds = tract_mask_vals==1\n",
    "        fa_tract = fa_vals[tract_mask_inds]\n",
    "        # Calculate average tract FA\n",
    "        mean_fa_tract = np.nanmean(fa_tract)\n",
    "        tract_vals[tract].append(mean_fa_tract)\n",
    "        \n",
    "        # METHOD 2\n",
    "        # Calculate average tractometry FA\n",
    "        tract_vals2[tract].append(np.mean(fa_tractometry_vals[tract]))\n",
    "        \n",
    "np.save('tract_vals.npy',tract_vals)\n",
    "np.save('tract_vals2.npy',tract_vals2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0bda5dd",
   "metadata": {},
   "source": [
    "# 7) Prepare pandas dataframe for statistical analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c60e136",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame()\n",
    "df['subjects'] = subs\n",
    "df['TR_RD'] = [int(i) for i in inds_rd]\n",
    "df['HT_LT'] = [int(i) for i in inds_low_t] \n",
    "df['RD_TR'] = [int(not(i)) for i in inds_rd] # Create reverted group labels to make violinplots look better\n",
    "df['LT_HT'] = [int(not(i)) for i in inds_low_t] # Create reverted group labels to make violinplots look better\n",
    "df['AGE'] = age\n",
    "df['SEX'] = sex\n",
    "df['HTLTxSEX'] = df['HT_LT']*df['SEX']\n",
    "df['TRRDxSEX'] = df['TR_RD']*df['SEX']\n",
    "df['HAND'] = hand\n",
    "df['ICV'] = icv\n",
    "df['TBV'] = tbv\n",
    "df['gFA'] = gFA\n",
    "df['TOWRE'] = towre\n",
    "df['SWE'] = swe\n",
    "df['PDE'] = pde\n",
    "df['all'] = 80 # Hack to make the violinplot work\n",
    "\n",
    "for tract in tracts:\n",
    "    df[tract] = tract_vals[tract]\n",
    "    df[tract+'2'] = tract_vals2[tract]\n",
    "    \n",
    "df.to_pickle('df.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62ec7fc4-26a5-485d-8266-a241c24ae81e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale = 2)\n",
    "sns.set_style(\"white\") # Format background style\n",
    "# Get dataframe with just covariates\n",
    "df_covars = df[['AGE','SEX','HAND','gFA','ICV','TBV','TOWRE']]\n",
    "\n",
    "# Calculate correlations across all pairwise columns\n",
    "corr = df_covars.corr(method='spearman')\n",
    "\n",
    "# Generate a mask for the upper triangle\n",
    "mask = np.triu(np.ones_like(corr, dtype=bool))\n",
    "\n",
    "# Set up the matplotlib figure\n",
    "f, ax = plt.subplots(figsize=(11, 9))\n",
    "\n",
    "# Generate a custom diverging colormap\n",
    "cmap = sns.diverging_palette(230, 20, as_cmap=True)\n",
    "\n",
    "# Draw the heatmap with the mask and correct aspect ratio\n",
    "sns.heatmap(corr, mask=mask, cmap=cmap, vmax = .45, vmin = -.45, center=0,\n",
    "            square=True, linewidths=.5, cbar_kws={\"shrink\": .5}, annot=True, annot_kws={\"fontsize\":15})\n",
    "#plt.savefig('covars_corr.eps',format='eps')\n",
    "plt.show()\n",
    "pingouin.rcorr(df_covars,method='spearman',padjust='fdr_bh')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nipype",
   "language": "python",
   "name": "nipype"
  },
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}