0.1_dataPrep_BLS-38.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adult brain measure estimation from BLS data - BLS-38\n",
    "As BLS is not publicly available due to ethical considerations, no data is included in this repository. The code is provided as a reference for the analysis of the BLS data. \n",
    "\n",
    "Since not all participants had a longitudinal scan yet, BLS-26 and BLS-38 were preprocessed in two different repositories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data preparation\n",
    "This whole notebook can be run for cortical thickness (CT) or surface area (SA). Both measures were used in the original publication. Adjust the variable `brain_measure` to `CT` or `SA` accordingly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from os.path import join\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import ttest_ind\n",
    "\n",
    "from abagen import fetch_desikan_killiany\n",
    "from neuroCombat import neuroCombat\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "# import custom code\n",
    "import sys\n",
    "sys.path.append('code')\n",
    "from preprocessing import load_freesurfer_aparc, load_freesurfer_aseg, calculate_scanner_difference, reorder_vars\n",
    "from euler import euler_hemi_combination, get_euler_outliers\n",
    "\n",
    "bls_data_dir = 'data/BLS'\n",
    "bls_freesurfer_outputs = join(bls_data_dir, 'freesurfer')  # where all freesurfer outputs are stored\n",
    "qc_dir_bls38 = join(bls_data_dir, 'euler_BLS-38')\n",
    "bls_out = join(bls_data_dir, 'derivatives')  # where the adapted data will be stored\n",
    "os.makedirs(bls_out, exist_ok=True)\n",
    "\n",
    "\n",
    "# adjust brain measurement for which the code should be run\n",
    "brain_measure = 'CT'  # CT or SA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load MRI related data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "110 participants in BLS had a longitundinal scan already and are included in the analysis\n"
     ]
    }
   ],
   "source": [
    "# brain_measure\n",
    "if brain_measure == 'CT':\n",
    "    mri_data_lh = load_freesurfer_aparc(file = join(bls_freesurfer_outputs, 'BLS-38_thickness_fs7.3.2_lh.txt'))\n",
    "    mri_data_rh = load_freesurfer_aparc(file = join(bls_freesurfer_outputs, 'BLS-38_thickness_fs7.3.2_rh.txt'))\n",
    "    mri_data = mri_data_lh.merge(mri_data_rh, on='participant')\n",
    "    #mri_data.drop(columns=['lh_MeanThickness_thickness', 'rh_MeanThickness_thickness'], inplace=True)\n",
    "    mri_data.rename(columns={'lh_MeanThickness_thickness':'lh_mean_CT2', 'rh_MeanThickness_thickness':'rh_mean_CT2'}, inplace=True)\n",
    "    bls_idps_idx = mri_data.filter(regex='_thickness').columns.to_list()\n",
    "    \n",
    "elif brain_measure == 'SA':\n",
    "    mri_data_lh = load_freesurfer_aparc(file = join(bls_freesurfer_outputs, 'BLS-38_surfaceArea_fs7.3.2_lh.txt'))\n",
    "    mri_data_rh = load_freesurfer_aparc(file = join(bls_freesurfer_outputs, 'BLS-38_surfaceArea_fs7.3.2_rh.txt'))\n",
    "    mri_data = mri_data_lh.merge(mri_data_rh, on='participant')\n",
    "    mri_data.drop(columns=['lh_WhiteSurfArea_area', 'rh_WhiteSurfArea_area'], inplace=True)\n",
    "    bls_idps_idx = mri_data.filter(regex='_area').columns.to_list()\n",
    "    \n",
    "else:\n",
    "    raise ValueError('brain_measure has to be either CT or SA')\n",
    "\n",
    "\n",
    "# aseg\n",
    "aseg = load_freesurfer_aseg(file = join(bls_freesurfer_outputs, 'BLS-38_aseg_fs7.3.2.txt'))\n",
    "aseg = aseg[['participant', 'EstimatedTotalIntraCranialVol','TotalGrayVol','SubCortGrayVol','CerebralWhiteMatterVol']]    \n",
    "\n",
    "# meta\n",
    "meta = pd.read_csv(join(bls_data_dir, 'BLS-38_meta_cleaned.csv'))\n",
    "meta = meta[['participant', 'blsgroup','phenotype','age', 'sex','sex_code']]\n",
    "\n",
    "# add additional information from BLS-26 meta\n",
    "meta_26 = pd.read_csv(join(bls_data_dir, 'BLS-26_meta_cleaned.csv'))\n",
    "meta_26 = meta_26[['participant', 'GA', 'BW']]\n",
    "\n",
    "# merge\n",
    "bls_38 = meta.merge(meta_26, on='participant')\n",
    "bls_38 = bls_38.merge(mri_data, on='participant')\n",
    "bls_38 = bls_38.merge(aseg, on='participant')\n",
    "\n",
    "# sort data\n",
    "bls_38.sort_values(by='participant', inplace=True)\n",
    "print(f'{bls_38.shape[0]} participants in BLS had a longitundinal scan already and are included in the analysis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of IDPs: 68\n"
     ]
    }
   ],
   "source": [
    "# access DK atlas and save idp_labels\n",
    "atlas = fetch_desikan_killiany(surface=True)\n",
    "atlas = pd.read_csv(atlas['info'])\n",
    "atlas = atlas[(atlas['structure'] == 'cortex') & (atlas['hemisphere'] == 'L')]\n",
    "atlas_labels_l = ['L_'+label for label in atlas['label']]\n",
    "atlas_labels_r = ['R_'+label for label in atlas['label']]\n",
    "\n",
    "desikan_idps = atlas_labels_l + atlas_labels_r\n",
    "print(f\"Number of IDPs: {len(bls_idps_idx)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename columns\n",
    "bls_38 = bls_38.rename(columns=dict(\n",
    "        SubCortGrayVol=\"sGMV\",\n",
    "        TotalGrayVol=\"GMV\",\n",
    "        EstimatedTotalIntraCranialVol=\"eTIV\",\n",
    "        CerebralWhiteMatterVol=\"WMV\"\n",
    "))\n",
    "bls_38 = bls_38.rename(columns=dict(zip(bls_idps_idx, desikan_idps)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality control using Euler Index\n",
    "Freesurfer's Euler number is a measure of the topological complexity of the reconstructed cortical surface. The Euler number is expected to be 2 for a closed surface. Freesurfer's Euler number has previously been used to identify subjects with poor quality reconstructions. Some informative publications are listed below:\n",
    "- [Rosen 2018](https://doi.org/10.1016/j.neuroimage.2017.12.059) \n",
    "- [Bethlehem 2022](https://doi.org/10.1038/s41586-022-04554-y)\n",
    "\n",
    "The Euler Index can be extracted with *mris_euler_index* from the recon-all output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median of Euler indices: -18.0\n",
      "Median absolute deviation (MAD): 8.0\n",
      "Euler index histogram saved to: data/BLS/euler_BLS-38/euler_histogram.png\n"
     ]
    }
   ],
   "source": [
    "# combine Euler Indices across hemispheres\n",
    "euler_hemi_combination(euler_file=os.path.join(qc_dir_bls38, 'QC_output_BLS-38.txt'), outdir=qc_dir_bls38, dataset_name='BLS-38')\n",
    "\n",
    "# calculate outliers based on Euler Indices\n",
    "outliers = get_euler_outliers(os.path.join(qc_dir_bls38, 'QC_output_combined_BLS-38.txt'), outdir=qc_dir_bls38, threshold=4.5)\n",
    "\n",
    "# adapt excluded subjects after manual QC\n",
    "outliers.append('BEST-BN-009')  # this subject does not have information for gestational age and was thus excluded\n",
    "\n",
    "# get outlier complete df\n",
    "outlier_data = bls_38[bls_38['participant'].isin(outliers)]\n",
    "outlier_data['blsgroup']\n",
    "\n",
    "with open(join(bls_data_dir, f'BLS-38_{brain_measure}_excluded_subjects.txt'), 'w') as f:\n",
    "    for item in outliers:\n",
    "        f.write(\"%s\\n\" % item)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data shape before exclusion: (110, 82)\n",
      "New size of data: (105, 82)\n"
     ]
    }
   ],
   "source": [
    "# actually exclude the outliers\n",
    "print(\"Data shape before exclusion:\", bls_38.shape)\n",
    "bls_38 = bls_38[~bls_38['participant'].isin(outliers)]\n",
    "\n",
    "print(\"New size of data:\", bls_38.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save data before hemi averaging\n",
    "# bls_38.to_csv(join(bls_out, f'BLS-38_{brain_measure}_raw.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average hemispheres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "regions = [r[2:] for r in desikan_idps[:34]]\n",
    "\n",
    "# Combine L_ and R_ values\n",
    "for region in regions:\n",
    "    bls_38[f'{brain_measure}_{region}'] = bls_38[[f'L_{region}', f'R_{region}']].mean(axis=1)\n",
    "    # drop L_ and R_ columns\n",
    "    bls_38 = bls_38.drop(columns=[f'L_{region}', f'R_{region}'])\n",
    "\n",
    "# create region list bilateral    \n",
    "desikan_idps_bilateral = [f'{brain_measure}_{region}' for region in regions]\n",
    "ctv_columns = ['GMV', 'WMV', 'sGMV', 'eTIV']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure data is ordered correctly\n",
    "bls_38 = bls_38.sort_values(by='participant')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adapt df for BrainChart framework\n",
    "BrainChart needs a certain format of the data. We will adapt the data accordingly. More information can be found [here](https://brainchart.shinyapps.io/brainchart/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# age\n",
    "bls_38.rename(columns={'age': 'Age'}, inplace=True)\n",
    "bls_38['age_days'] = (bls_38.Age * 365.245) + 280\n",
    "\n",
    "# sex\n",
    "bls_38['sex'] = bls_38['sex'].map({'female':'Female', 'male':'Male'})\n",
    "\n",
    "# dx\n",
    "bls_38['dx'] = bls_38['phenotype'].map({'term':'CN', 'preterm':'preterm'})\n",
    "\n",
    "# other columns\n",
    "bls_38['study'] = 'BLS-38'\n",
    "bls_38['fs_version'] = 'Custom'\n",
    "bls_38['country'] = 'Multisite'\n",
    "bls_38['run'] = 1\n",
    "bls_38['session'] = 1\n",
    "\n",
    "# drop columns\n",
    "bls_38.drop(columns=['blsgroup', 'phenotype', 'sex_code'], inplace=True)\n",
    "\n",
    "# reshape and rename last things\n",
    "bls_38_final = reorder_vars(['participant', 'Age', 'age_days', 'sex', 'study', 'fs_version','country', 'run', \n",
    "                            'session', 'dx'], bls_38, desikan_idps_bilateral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save\n",
    "bls_38_final.to_csv(join(bls_out, f'BLS-38_{brain_measure}_preprocessed.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stats\n",
    "Summary stats shown in Supp Table S1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overall number of unique subjects with longitudinal data: n = 105\n",
      "Preterm stats: n =  52\n",
      "Male      30\n",
      "Female    22\n",
      "Name: sex, dtype: int64\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Age</th>\n",
       "      <th>GA</th>\n",
       "      <th>BW</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>52.00</td>\n",
       "      <td>52.00</td>\n",
       "      <td>52.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>38.10</td>\n",
       "      <td>30.10</td>\n",
       "      <td>1318.56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.46</td>\n",
       "      <td>1.89</td>\n",
       "      <td>335.99</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>37.36</td>\n",
       "      <td>25.00</td>\n",
       "      <td>730.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>37.74</td>\n",
       "      <td>29.00</td>\n",
       "      <td>990.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>38.19</td>\n",
       "      <td>30.00</td>\n",
       "      <td>1345.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>38.41</td>\n",
       "      <td>31.00</td>\n",
       "      <td>1512.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>39.07</td>\n",
       "      <td>35.00</td>\n",
       "      <td>2070.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Age     GA       BW\n",
       "count  52.00  52.00    52.00\n",
       "mean   38.10  30.10  1318.56\n",
       "std     0.46   1.89   335.99\n",
       "min    37.36  25.00   730.00\n",
       "25%    37.74  29.00   990.00\n",
       "50%    38.19  30.00  1345.00\n",
       "75%    38.41  31.00  1512.50\n",
       "max    39.07  35.00  2070.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Overall number of unique subjects with longitudinal data: n =\", len(bls_38_final))\n",
    "\n",
    "bls_38_final_pt = bls_38_final[bls_38_final['dx'] == 'preterm']\n",
    "bls_38_final_cn = bls_38_final[bls_38_final['dx'] == 'CN']\n",
    "\n",
    "print('Preterm stats: n = ', len(bls_38_final_pt))\n",
    "print(bls_38_final_pt['sex'].value_counts())\n",
    "display(bls_38_final_pt[['Age', 'GA', 'BW']].describe().round(2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Full-term stats: n =  53\n",
      "Male      30\n",
      "Female    23\n",
      "Name: sex, dtype: int64\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Age</th>\n",
       "      <th>GA</th>\n",
       "      <th>BW</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>53.00</td>\n",
       "      <td>53.00</td>\n",
       "      <td>53.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>38.21</td>\n",
       "      <td>39.89</td>\n",
       "      <td>3474.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.40</td>\n",
       "      <td>1.09</td>\n",
       "      <td>419.71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>37.41</td>\n",
       "      <td>37.00</td>\n",
       "      <td>2600.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>37.93</td>\n",
       "      <td>39.00</td>\n",
       "      <td>3200.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>38.22</td>\n",
       "      <td>40.00</td>\n",
       "      <td>3450.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>38.55</td>\n",
       "      <td>41.00</td>\n",
       "      <td>3730.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>38.91</td>\n",
       "      <td>42.00</td>\n",
       "      <td>4670.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Age     GA       BW\n",
       "count  53.00  53.00    53.00\n",
       "mean   38.21  39.89  3474.89\n",
       "std     0.40   1.09   419.71\n",
       "min    37.41  37.00  2600.00\n",
       "25%    37.93  39.00  3200.00\n",
       "50%    38.22  40.00  3450.00\n",
       "75%    38.55  41.00  3730.00\n",
       "max    38.91  42.00  4670.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Full-term stats: n = ', len(bls_38_final_cn))\n",
    "print(bls_38_final_cn['sex'].value_counts())\n",
    "display(bls_38_final_cn[['Age', 'GA', 'BW']].describe().round(2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}