univariate_vlsm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Univariate VLSM\n",
    "\n",
    "This script can be used to perform univariate voxel-based lesion-symptom mapping (VLSM) or be used as feature-selection for multi-variate LSM (e.g. with support-vector regression (SVR)).\n",
    "\n",
    "There are two assumptions on your data:\n",
    "1. You have binary lesions maps (0=background, 1=lesion) in MNI space\n",
    "2. The output is a continuous variable (e.g. Z-scores)\n",
    "\n",
    "The script performs the following steps:\n",
    "1. A lesion prevalence analysis, showing how many subjects have a lesion in a specific voxel. This is used to only assess voxels where a sufficient number of subjects have a lesion (e.g. >= 5 subjects)\n",
    "2. A voxel-wise t-test\n",
    "3. Computing a power-map per voxel\n",
    "  1. Calculating the effect-size per voxel and taking the 99th percentile as the overall (fixed) effect size\n",
    "  2. Computing the power-map per voxel, using the fixed effect size from (a)\n",
    "4. Performing multiple testing correction\n",
    "  1. FDR correction (Benjamini/Hochberg)\n",
    "  2. Optionally: a permutation minT/maxP correction\n",
    "5. Saving everything as nifti images in MNI space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import json\n",
    "import os\n",
    "\n",
    "from joblib import Parallel, delayed\n",
    "import numpy as np\n",
    "import pandas\n",
    "from scipy.stats import ttest_ind\n",
    "import SimpleITK as sitk\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "from statsmodels.stats.power import tt_ind_solve_power\n",
    "from tqdm.notebook import trange, tqdm\n",
    "from tqdm.contrib import tenumerate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters\n",
    "\n",
    "All relevant parameters for this script can be set in the box below, related to either the input, settings, or output.\n",
    "\n",
    "### Input \n",
    "\n",
    "- **dir_lesion:** the directory that contains the binary lesion maps per subject\n",
    "- **design_document:** spreadsheet containing two or more columns: *first column:* name of the lesion map for every subject (should exist in dir_lesion or a subfolder (see below)); *second column:* continuous output variable, e.g. a Z-score on a specific domain; *optionally* more columns with additional output variables\n",
    "- **data_in_subfolders:** specifies whether the data is directly in *dir_lesion* (False) or whether it is in subfolders (True). In the latter case, the filename in the first column of the *design_document* should be: COHORT_SUBJECT.nii.gz and this file is located in *dir_lesion*/COHORT/COHORT_SUBJECT.nii.gz  \n",
    "- **domain:** which column in the design_document to use as output (note: 0-based index). By default: 1.\n",
    "\n",
    "### Settings\n",
    "\n",
    "- **subject_threshold:** minimum number of subjects with a lesion in a voxel\n",
    "- **alpha:** alpha value to use in t-test and multiple testing correction\n",
    "- **perform_multiple_testing_permutation:** whether to perform multiple testing correction using permutations (minT/maxP ; True) or not (False). Note: this is extremely slow and can easily take more than 24 hours.\n",
    "- **num_permutations:** number of permutations for the multiple testing correction\n",
    "- **n_jobs:** number of parallel computation jobs (set at 1 when unsure)\n",
    "\n",
    "### Output\n",
    "- **output_base_path:** output directory into which all results will be written\n",
    "- **output_name_\\*:** output filenames for the various maps that are computed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input data\n",
    "dir_lesion = r\"\"\n",
    "design_document = r\"\"\n",
    "data_in_subfolders = True\n",
    "domain = 1\n",
    "\n",
    "# Settings\n",
    "subject_threshold = 5\n",
    "alpha = 0.05\n",
    "perform_multiple_testing_permutation = False  # Note: super slow !\n",
    "num_permutations = 1000\n",
    "n_jobs = 35\n",
    "\n",
    "# Output\n",
    "output_base_path = r\"\"\n",
    "\n",
    "# Output for lesion prevalence map\n",
    "output_name_lesion_prevalence = output_base_path + \"_lesion_prevalence.nii\"\n",
    "\n",
    "# Output for t-test and multiple testing correction\n",
    "output_name_tdata = output_base_path + \"_tmap.nii\"\n",
    "output_name_pdata = output_base_path + \"_1-p_not_corrected.nii\"\n",
    "output_name_pcorrdata = output_base_path + \"_1-p_fdr_corrected.nii\"\n",
    "\n",
    "# Output for power calculation\n",
    "output_name_powermap = output_base_path + \"_power_map.nii\"\n",
    "output_name_effectmap = output_base_path + \"_effect_map.nii\"\n",
    "\n",
    "# Output for permutation tests (optional, depending on perform_multiple_testing_permutation==True)\n",
    "output_name_ppermcorrdata = output_base_path + \"_1-p_permutation_corrected.nii\"\n",
    "output_name_tstatistics = output_base_path + \"_tstatistics.dat\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do not make changes below here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the design document and load the Z-scores\n",
    "df = pandas.read_excel(design_document, header=None)\n",
    "z_score = df.iloc[:, domain].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize lesion matrix and load all lesion data\n",
    "lesion_filename = df[0][0]\n",
    "if data_in_subfolders:\n",
    "    subfolder = lesion_filename.split(\"_\")[0]\n",
    "    lesion_filename = os.path.join(subfolder, lesion_filename)\n",
    "nii = sitk.ReadImage(os.path.join(dir_lesion, lesion_filename))\n",
    "\n",
    "# The raw_lesion_matrix has the shape: number of subjects, number of voxels\n",
    "raw_lesion_matrix = np.zeros((len(df.index), sitk.GetArrayViewFromImage(nii).size), np.int8)\n",
    "for i, lesion_filename in tenumerate(df[0]):\n",
    "    if data_in_subfolders:\n",
    "        subfolder = lesion_filename.split(\"_\")[0]\n",
    "        lesion_filename = os.path.join(subfolder, lesion_filename)\n",
    "\n",
    "    nii = sitk.ReadImage(os.path.join(dir_lesion, lesion_filename))\n",
    "    raw_lesion_matrix[i] = sitk.GetArrayViewFromImage(nii).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the lesion prevalence map\n",
    "lesion_prevalence = np.sum(raw_lesion_matrix, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine which voxels to test (sufficient number of subjects per voxel)\n",
    "index_test_mask = lesion_prevalence >= subject_threshold\n",
    "index_test_voxels = np.argwhere(index_test_mask)\n",
    "n_test_voxels = len(index_test_voxels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform t-test\n",
    "t_data_v = np.zeros(n_test_voxels)\n",
    "pvalue_v = np.zeros(n_test_voxels)\n",
    "\n",
    "\n",
    "def do_ttest(i):\n",
    "    group_nc = raw_lesion_matrix[:, index_test_voxels[i][0]] == 0\n",
    "    group_damage = raw_lesion_matrix[:, index_test_voxels[i][0]] > 0\n",
    "    \n",
    "    return ttest_ind(z_score[group_nc], z_score[group_damage], equal_var=True)\n",
    "\n",
    "\n",
    "result = Parallel(n_jobs=n_jobs)(delayed(do_ttest)(i) for i in trange(n_test_voxels))\n",
    "t_data_v, pvalue_v = zip(*result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform multiple testing correction using permutations. Results are stored (per permutation) in\n",
    "# a text file. This way, you can resume and append the results.\n",
    "if perform_multiple_testing_permutation:\n",
    "    rng = np.random.default_rng()\n",
    "\n",
    "    for _ in trange(num_permutations):\n",
    "        permuted_z_score = rng.permutation(z_score)\n",
    "        \n",
    "        \n",
    "        def do_permuted_ttest(i):\n",
    "            group_nc = raw_lesion_matrix[:, index_test_voxels[i][0]] == 0\n",
    "            group_damage = raw_lesion_matrix[:, index_test_voxels[i][0]] > 0\n",
    "\n",
    "            return ttest_ind(permuted_z_score[group_nc], permuted_z_score[group_damage], equal_var=True)\n",
    "        \n",
    "        \n",
    "        permuted_result = Parallel(n_jobs=n_jobs)(delayed(do_permuted_ttest)(i) for i in trange(n_test_voxels, leave=False))\n",
    "        permuted_t_data_v, _ = zip(*permuted_result)\n",
    "\n",
    "        with open(output_name_tstatistics, 'ab') as tstatistics_file:\n",
    "            np.savetxt(tstatistics_file, np.array([max(permuted_t_data_v)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in the results from the previous cell and compute the multiple testing correction\n",
    "if perform_multiple_testing_permutation:\n",
    "    with open(output_name_tstatistics, 'rb') as tstatistics_file:\n",
    "        max_t_distribution = np.loadtxt(tstatistics_file)\n",
    "\n",
    "    pvalue_permutation_corrected = np.mean(max_t_distribution[:, np.newaxis] > t_data_v, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute all effect sizes\n",
    "z_score_std = np.std(z_score)\n",
    "def do_effect_size(i):\n",
    "    group_nc = raw_lesion_matrix[:, index_test_voxels[i][0]] == 0\n",
    "    group_damage = raw_lesion_matrix[:, index_test_voxels[i][0]] > 0    \n",
    "    \n",
    "    z_score_nc = np.mean(z_score[group_nc])\n",
    "    z_score_damage = np.mean(z_score[group_damage])\n",
    "    \n",
    "    return (z_score_nc - z_score_damage) / z_score_std\n",
    "\n",
    "result_effect = Parallel(n_jobs=n_jobs)(delayed(do_effect_size)(i) for i in trange(n_test_voxels))\n",
    "\n",
    "# We use the 99th percentile as a fixed effect size\n",
    "fixed_effect_size = np.percentile(result_effect, 99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute power maps. There are two functions:\n",
    "# do_tt_ind_solve_power: determine the effect size per voxel and compute the power\n",
    "# do_tt_ind_solve_power_fixed_effect_size: use the fixed effect size (computed above) for all voxels\n",
    "\n",
    "z_score_std = np.std(z_score)\n",
    "def do_tt_ind_solve_power(i):\n",
    "    group_nc = raw_lesion_matrix[:, index_test_voxels[i][0]] == 0\n",
    "    group_damage = raw_lesion_matrix[:, index_test_voxels[i][0]] > 0    \n",
    "    \n",
    "    z_score_nc = np.mean(z_score[group_nc])\n",
    "    z_score_damage = np.mean(z_score[group_damage])\n",
    "    effect_size = (z_score_nc - z_score_damage) / z_score_std\n",
    "\n",
    "    nobs1 = np.sum(group_nc)\n",
    "    ratio = np.sum(group_damage) / nobs1\n",
    "    \n",
    "    return tt_ind_solve_power(effect_size=effect_size, nobs1=nobs1, alpha=alpha, power=None, ratio=ratio, alternative='two-sided')\n",
    "\n",
    "\n",
    "def do_tt_ind_solve_power_fixed_effect_size(i):   \n",
    "    group_nc = raw_lesion_matrix[:, index_test_voxels[i][0]] == 0\n",
    "    group_damage = raw_lesion_matrix[:, index_test_voxels[i][0]] > 0    \n",
    "     \n",
    "    nobs1 = np.sum(group_nc)\n",
    "    ratio = np.sum(group_damage) / nobs1\n",
    "\n",
    "    return tt_ind_solve_power(effect_size=fixed_effect_size, nobs1=nobs1, alpha=alpha, power=None, ratio=ratio, alternative='two-sided')\n",
    "\n",
    "\n",
    "# Change the function in delayed() to the desired way of computing the power.\n",
    "result_power = Parallel(n_jobs=n_jobs)(delayed(do_tt_ind_solve_power_fixed_effect_size)(i) for i in trange(n_test_voxels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correction for multiple testing\n",
    "_, pvals_corrected, _, _ = multipletests(pvalue_v, alpha=alpha, method='fdr_bh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the MNI template, to use as a reference image for all output images\n",
    "ref_nii = sitk.ReadImage(os.path.join(os.getcwd(), 'src', 'Atlas', 'LSM_reference_1mm_MNI152.nii'))\n",
    "ref_nii_arr = sitk.GetArrayViewFromImage(ref_nii)\n",
    "\n",
    "# Save all maps as nii\n",
    "nii = sitk.GetImageFromArray(lesion_prevalence.reshape(ref_nii_arr.shape).astype(float))\n",
    "nii.CopyInformation(ref_nii)\n",
    "sitk.WriteImage(nii, output_name_lesion_prevalence)\n",
    "\n",
    "\n",
    "def save_result_image(data, output_name):\n",
    "    data_arr = np.zeros(ref_nii_arr.shape, np.float32)\n",
    "    data_arr[index_test_mask.reshape(ref_nii_arr.shape)] = data\n",
    "    data_img = sitk.GetImageFromArray(data_arr)\n",
    "    data_img.CopyInformation(ref_nii)\n",
    "    sitk.WriteImage(data_img, output_name)\n",
    "\n",
    "\n",
    "save_result_image(np.array(t_data_v), output_name_tdata)\n",
    "save_result_image(1-np.array(pvalue_v), output_name_pdata)\n",
    "save_result_image(1-pvals_corrected, output_name_pcorrdata)\n",
    "if perform_multiple_testing_permutation:\n",
    "    save_result_image(1-pvalue_permutation_corrected, output_name_ppermcorrdata)\n",
    "save_result_image(np.array(result_power), output_name_powermap)\n",
    "save_result_image(np.array(result_effect), output_name_effectmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: save description of input, settings, outputs, date/time, etc"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}