sEEG MRI co-localisation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This script requires the following for each participant\n",
    "# 1) Classifier output from the structural MRI analysis - this should be a probability map between 0 and 1. \n",
    "#(see FCDdetection scripts for how to run a neural network to detect FCDs)\n",
    "# 2) A .fcsv file for each electrode implanted. This should contain the co-ordinates of each contact on the electrode\n",
    "# 3) A spreadsheet of participants, electrodes implanted, ictal contacts (see example on github)\n",
    "\n",
    "#This script does the following:\n",
    "#1. Load in the classifier output from the classifier training cohort and find the optimal threshold to threshold the probability map\n",
    "#2. Calculates the sensitivity of the classifier on the training cohort at this threshold\n",
    "#3. Load in the classifier output tested on the controls and finds out the specificity at this optimal threshold\n",
    "#4. Load in the SEEG patients' classifier output, threshold and cluster their outputs \n",
    "#5. Save an interactive plot as a .html file of the electrodes and the clusters in the SEEG patients\n",
    "#6. Plot the relationship of the ictal electrodes to the top cluster\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary packages\n",
    "import io_meld as io\n",
    "import numpy as np\n",
    "import os\n",
    "import mesh_tools as mt\n",
    "from matplotlib import pyplot as plt\n",
    "import glob\n",
    "%matplotlib inline\n",
    "import subprocess\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "from nilearn import plotting, surface,datasets\n",
    "%matplotlib notebook\n",
    "from matplotlib import cm as mpl_cm\n",
    "from nilearn.plotting import cm\n",
    "import matplotlib as mpl\n",
    "from matplotlib.colors import Normalize, LinearSegmentedColormap\n",
    "import io_mesh\n",
    "import json\n",
    "from scipy import sparse\n",
    "from scipy import ndimage\n",
    "from nilearn import datasets\n",
    "import nilearn\n",
    "import numbers\n",
    "from nilearn.plotting.js_plotting_utils import (add_js_lib, HTMLDocument, mesh_to_plotly,\n",
    "                                encode, colorscale, get_html_template,\n",
    "                                to_color_strings)\n",
    "import argparse\n",
    "\n",
    "from nilearn.plotting.js_plotting_utils import (add_js_lib, HTMLDocument, mesh_to_plotly,\n",
    "                                encode, colorscale, get_html_template,\n",
    "                                to_color_strings)\n",
    "from nilearn.surface import load_surf_data, load_surf_mesh, vol_to_surf\n",
    "\n",
    "from nilearn._utils.compat import _basestring\n",
    "import PtitPrince as pt\n",
    "import seaborn as sns "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define functions used\n",
    "def get_neighbours_from_tris(tris, label=None):\n",
    "    \"\"\"Get surface neighbours from tris\n",
    "        Input: tris\n",
    "         Returns Nested list. Each list corresponds \n",
    "        to the ordered neighbours for the given vertex\"\"\"\n",
    "    n_vert=np.max(tris)+1\n",
    "    neighbours=[[] for i in range(n_vert)]\n",
    "    for tri in tris:\n",
    "        neighbours[tri[0]].extend([tri[1],tri[2]])\n",
    "        neighbours[tri[2]].extend([tri[0],tri[1]])\n",
    "        neighbours[tri[1]].extend([tri[2],tri[0]])\n",
    "    #Get unique neighbours\n",
    "    for k in range(len(neighbours)):      \n",
    "        if label is not None:\n",
    "            neighbours[k] = set(neighbours[k]).intersection(label)\n",
    "        else :\n",
    "            neighbours[k]=f7(neighbours[k])\n",
    "    return neighbours;\n",
    "\n",
    "\n",
    "def f7(seq):\n",
    "    \"\"\"returns uniques but in order of original appearance.\n",
    "    Used to retain neighbour triangle relationship\"\"\"\n",
    "    seen = set()\n",
    "    seen_add = seen.add\n",
    "    return [x for x in seq if not (x in seen or seen_add(x))];\n",
    "\n",
    "def flatten(nested_list):\n",
    "        \"\"\"flatten nested list\"\"\"\n",
    "        #test if nested\n",
    "        if any(isinstance(i, list) for i in nested_list):\n",
    "            flat_list = [item for sublist in nested_list for item in sublist]\n",
    "            return flat_list\n",
    "        else:\n",
    "            return nested_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define functions used\n",
    "def cluster_above_threshold(overlay, neighbours,threshold):\n",
    "    \"\"\"cluster outputs above threshold and returns cluster index at vertex locations\"\"\"\n",
    "    unranked=overlay>=threshold\n",
    "    clust_ind=0\n",
    "    clusters=np.zeros(len(overlay))\n",
    "    while np.sum(unranked)>0:\n",
    "        cluster=np.zeros(len(overlay))\n",
    "        clust_ind+=1\n",
    "        new_verts=np.argmax(unranked)\n",
    "        unranked[new_verts]=0\n",
    "        cluster[new_verts]=clust_ind\n",
    "        old_len=np.sum(unranked)+1\n",
    "        while np.sum(unranked)<old_len:\n",
    "            neighbours_vert=np.zeros(len(overlay))\n",
    "            try :\n",
    "                neighbours_vert[flatten(neighbours[new_verts])]=1\n",
    "            except TypeError:\n",
    "                neighbours_vert[neighbours[new_verts]]=1\n",
    "            new_verts=np.logical_and(neighbours_vert,unranked)\n",
    "            old_len=np.sum(unranked)\n",
    "            cluster[new_verts]=clust_ind\n",
    "            unranked[new_verts]=0\n",
    "        clusters+=cluster\n",
    "    return clusters.astype(int)\n",
    "\n",
    "def filter_cluster(clusters,areas, minimum_area):\n",
    "    \"\"\"filter out clusters smaller than minimum areas\"\"\"\n",
    "    #get existing clusters\n",
    "    clusters_indices=np.unique(clusters)\n",
    "    #start after the zero index\n",
    "    filtered_clusters=clusters.copy()\n",
    "    #store which vertices grew\n",
    "    for c in clusters_indices[1:]:\n",
    "        #calculate cluster area\n",
    "        area=np.sum(areas[clusters==c])\n",
    "        #if cluster area below threshold, set all vertices to zero\n",
    "        if area < minimum_area:\n",
    "            filtered_clusters[clusters==c]=0\n",
    "    return filtered_clusters\n",
    "\n",
    "def clump_clusters(clusters,neighbours, n_steps=2):\n",
    "    \"\"\"group clusters separated by n_steps\"\"\"\n",
    "    cluster_indices=np.unique(clusters)\n",
    "    grouped_clusters=clusters.copy()\n",
    "    growth_rings=np.zeros((len(cluster_indices),len(clusters)))\n",
    "    for c in cluster_indices[1:]:\n",
    "        cluster=clusters==c\n",
    "        #grow out cluster\n",
    "        for step in range(n_steps):\n",
    "            cluster=np.unique(flatten(neighbours[cluster]))\n",
    "        #check if grown cluster contains new indices from nearby clusters\n",
    "        indices_in_grown_cluster=np.unique(clusters[cluster])\n",
    "        if len(indices_in_grown_cluster)>2:\n",
    "            #set those clusters to that value\n",
    "            growth_rings[c,cluster]=1\n",
    "            for index in indices_in_grown_cluster[1:]:\n",
    "                #look up in grouped clusters to pass forward any changed values from previous iterations\n",
    "                if grouped_clusters[clusters==c][0] >0:\n",
    "                    new_value = grouped_clusters[clusters==c][0]\n",
    "                    sum_growth=growth_rings[new_value,:]+growth_rings[c,:]\n",
    "                    link=sum_growth==2\n",
    "                    #print(link)\n",
    "                else:\n",
    "                    new_value= c\n",
    "                grouped_clusters[clusters==index]=new_value\n",
    "               # grouped_clusters[link]=new_value\n",
    "    return grouped_clusters\n",
    "            #save growth ring\n",
    "            \n",
    "\n",
    "from scipy.stats import skew\n",
    "def rank_clusters(clusters,overlay,ranking=\"mean\"):\n",
    "    \"\"\"reorder indexing to order cluster by mean/max value \"\"\"\n",
    "    clusters_indices=np.unique(clusters)[1:]\n",
    "    values=np.zeros(len(clusters_indices))\n",
    "    for k,c in enumerate(clusters_indices):\n",
    "        #calculate cluster area\n",
    "        if ranking==\"mean\":\n",
    "            values[k]=np.mean(overlay[clusters==c])\n",
    "        elif ranking==\"max\":\n",
    "            values[k]=np.max(overlay[clusters==c])\n",
    "        elif ranking==\"skew\":\n",
    "            values[k]=skew(overlay[clusters==c])\n",
    "    ranked_clusters=clusters.copy()\n",
    "    order=np.argsort(values)\n",
    "    for k,i in enumerate(order):\n",
    "        ranked_clusters[clusters==clusters_indices[k]]=i+1\n",
    "    return ranked_clusters\n",
    "\n",
    "def check_cluster_overlap_with_lesion_mask(clusters, mask, cluster_index=1):\n",
    "    \"\"\"test if lesion mask overlaps with a specific cluster\n",
    "    cluster_index - default tests just 1 top cluster overlap\"\"\"\n",
    "    return np.logical_and(clusters==cluster_index,mask==1).any()\n",
    "\n",
    "def check_any_overlap_with_lesion_mask(clusters, mask):\n",
    "    \"\"\"test if lesion mask overlaps with any of the clusters\"\"\"\n",
    "    return np.logical_and(clusters>0,mask==1).any()\n",
    "\n",
    "\n",
    "def return_cluster_overlap_with_lesion_mask(clusters, mask):\n",
    "    \"\"\"test if lesion mask overlaps with any of the clusters\n",
    "    cluster_index - default tests just 1 top cluster overlap\"\"\"\n",
    "    overlaps=np.unique(clusters[np.logical_and(clusters>0,mask==1)])\n",
    "    if overlaps.any():\n",
    "        return overlaps\n",
    "    else: \n",
    "        return [0]\n",
    "    \n",
    "def get_optimal_threshold(sensitivity, specificity, thresholds, resolution=10):\n",
    "    \"\"\"calculate optimal sensitivity\"\"\"\n",
    "    y=sensitivity+specificity\n",
    "    x=thresholds\n",
    "    upsampled=np.repeat(y,10)\n",
    "    smoothed=ndimage.filters.gaussian_filter1d(upsampled,resolution)\n",
    "    upsampled_thresholds=np.linspace(0,1,len(thresholds)*resolution)\n",
    "    return np.round(upsampled_thresholds[np.argmax(smoothed)],2)\n",
    "\n",
    "\n",
    "def defrag_surface(clusters,neighbours,steps=5):\n",
    "    \"\"\" Removes small holes in clusters by expanding and shrinking the cluster\"\"\"\n",
    "    #find basic lesion vertics\n",
    "    cluster_indices=np.unique(clusters)[1:]\n",
    "    clusters_defragged=clusters.copy()\n",
    "    for cluster in cluster_indices:\n",
    "        patch=clusters==cluster\n",
    "        #expanding stage. add neighbours\n",
    "        for k in range(steps):\n",
    "            patch[flatten(neighbours[patch])]=True\n",
    "        #shrinking stage. remove edge vertices\n",
    "        not_patch=patch==False\n",
    "        for k in range(steps):\n",
    "            not_patch[flatten(neighbours[not_patch])]=True\n",
    "        clusters_defragged[not_patch==False]=cluster\n",
    "    return clusters_defragged;\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def label2object(outname,grey,white,label,shrink=False):\n",
    "    \"\"\"reduce object to just the label vertices\n",
    "    NB if there are isolated, unconnected vertices this causes bugs\n",
    "    shrink - removes non-label vertices\"\"\"\n",
    "    gsurf=io_mesh.load_mesh_geometry(grey)\n",
    "    wsurf=io_mesh.load_mesh_geometry(white)\n",
    "    label_ones=np.zeros(len(gsurf['coords']))\n",
    "    label_ones[label]=1\n",
    "    #extract only triangles inside label\n",
    "    labeltris=[]\n",
    "    for tri in gsurf['faces']:\n",
    "        if sum(label_ones[tri]) == 3:\n",
    "            labeltris.append(tri.tolist())\n",
    "    label_surf={}\n",
    "    if shrink:\n",
    "        labeltris=flatten(labeltris)\n",
    "        newlabeltris=[]\n",
    "        #create new indexing\n",
    "        labellist=sorted(label.tolist())\n",
    "        for l in labeltris:\n",
    "            newlabeltris.append(labellist.index(l))\n",
    "        newlabeltris=np.reshape(newlabeltris,(len(newlabeltris)/3,3))\n",
    "        #create new object\n",
    "        label_surf['coords']=np.vstack((gsurf['coords'][label],wsurf['coords'][label]))\n",
    "        label_surf['faces']=np.vstack((np.array(newlabeltris),np.array(newlabeltris)+len(label)))\n",
    "    else:\n",
    "        label_surf['coords']=np.vstack((gsurf['coords'],wsurf['coords']))\n",
    "        label_surf['faces']=np.vstack((np.array(labeltris),np.array(labeltris)+len(gsurf['coords'])))\n",
    "    io_mesh.save_mesh_geometry(outname,label_surf)\n",
    "    return\n",
    "\n",
    "def tidy_integer_labels(label,neighbours):\n",
    "    \"\"\" tidy overlay intended to be integers, but with interpolation errors\n",
    "    find all values greater than zero. if connected to higher value, set as higher value\"\"\"\n",
    "    vertices=np.where(label>0)[0]\n",
    "    change=1\n",
    "    #count number fo vertices that have changed, keep iterating until change stops\n",
    "    while change>0:\n",
    "        #store the sum of label values from previous iteration\n",
    "        old_sum=np.sum(label)\n",
    "        #for each non-zero vertex...\n",
    "        for vertex in vertices:\n",
    "            #find the highest value from among the neighbours\n",
    "            max_n=np.max(label[neighbours[vertex]])\n",
    "            #if the highest neighbouring value is above the value of the vertex, replace with neighbour value\n",
    "            if max_n>=label[vertex] :\n",
    "                label[vertex]=max_n\n",
    "            #otherwise set it to zero\n",
    "            else:\n",
    "                label[vertex]=0\n",
    "        change=np.sum(label)-old_sum\n",
    "        vertices=np.where(label>0)[0]\n",
    "    non_ints=np.where(label%1>0)[0]\n",
    "    label[non_ints]=0\n",
    "    return label.astype(int)\n",
    "        \n",
    "    \n",
    "def lesion_coordinates(grey,white,label_name,cluster_no=1):\n",
    "    \"\"\"returns coordinates of grey and white surfaces of lesion defined by label\"\"\"\n",
    "    gsurf=io_mesh.load_mesh_geometry(grey)\n",
    "    wsurf=io_mesh.load_mesh_geometry(white)\n",
    "    neighbours=get_neighbours_from_tris(gsurf['faces'])\n",
    "    label=tidy_integer_labels(np.array(io.import_mgh(label_name)),neighbours)    \n",
    "    lesion_coords=np.vstack((gsurf['coords'][label==cluster_no],wsurf['coords'][label==cluster_no]))\n",
    "    return lesion_coords\n",
    "\n",
    "# Define functions for seeg data\n",
    "def get_fiducials(subject_id):\n",
    "    \"\"\"get co-ordinates of electrode contacts 'fiducials' from a dataframe \"\"\"\n",
    "    fiducials=np.array(Fiducials[subject_id]).astype(str)\n",
    "    fiducials = fiducials[fiducials != 'nan']\n",
    "    fiducials = fiducials[fiducials != 'NaN']\n",
    "    return fiducials\n",
    "\n",
    "def get_electrodes(sheet, subject_id):\n",
    "    \"\"\" Retrieves the list of electrodes from the spreadsheet\"\"\"\n",
    "    row = np.array(sheet.loc[sheet['ID'] == subject_id]).astype(str)[0][1:]\n",
    "    row = row[row!='nan']\n",
    "    return row\n",
    "\n",
    "def closest_node(node, nodes):\n",
    "    deltas = nodes - node\n",
    "    dist_2 = np.einsum('ij,ij->i', deltas, deltas)\n",
    "    return np.argmin(dist_2), np.sqrt(np.min(dist_2))\n",
    "\n",
    "def _get_markers(coords, colors):\n",
    "    \"\"\" create connectome dictionary containing markers (electrode coordinates)\n",
    "    Colour each marker by colors, an array specifiying the desired colours.  \"\"\"\n",
    "    connectome = {}\n",
    "    coords = np.asarray(coords, dtype='<f4')\n",
    "    x, y, z = coords.T\n",
    "    for coord, cname in [(x, \"x\"), (y, \"y\"), (z, \"z\")]:\n",
    "        connectome[\"_con_{}\".format(cname)] = encode(\n",
    "            np.asarray(coord, dtype='<f4'))\n",
    "    connectome[\"marker_color\"] = to_color_strings(colors)\n",
    "    connectome[\"markers_only\"] = True\n",
    "    return connectome\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the minimum cluster area (originally 50mm2)\n",
    "min_area_threshold=50\n",
    "# Choose the default way to assess the top cluster - can be mean, max or skew\n",
    "clusters_by_mean_max=\"skew\"\n",
    "# Choose how many thresholds to test to pick optimal threshold - originally 21\n",
    "thresholds=np.linspace(0,1,21)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Path to training subjects directory\n",
    "pat_subjects_dir='/path/to/train_dir'\n",
    "#Path to control subjects directory\n",
    "cont_subjects_dir='/path/to/control_dir'\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates list of training subjects ids - looks for patient ids beginning \"FCD\"\n",
    "subject_ids=glob.glob(os.path.join(pat_subjects_dir,'FCD*'))\n",
    "# List of patients to remove from subject IDs\n",
    "Remove=['FCD_patient_1', 'FCD_patient_2']\n",
    "patient_ids=[]\n",
    "for subject in subject_ids:\n",
    "    split_sub=subject.split('/')[-1]\n",
    "    if split_sub not in Remove:\n",
    "        patient_ids.append(split_sub)\n",
    "        \n",
    "#Save the list of training IDs        \n",
    "np.savetxt('patient_ids_sEEG.txt',np.array(patient_ids), fmt = '%s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creates list of control ids - looks for IDs beginning C_\n",
    "subject_ids=glob.glob(os.path.join(cont_subjects_dir,'C_*'))\n",
    "control_ids=[]\n",
    "for subject in subject_ids:\n",
    "    split_sub=subject.split('/')[-1]\n",
    "    if split_sub not in Remove:\n",
    "        control_ids.append(split_sub)\n",
    "#Saves the list of control IDs as a .txt file\n",
    "np.savetxt('control_ids_sEEG.txt',np.array(control_ids), fmt = '%s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#preload fsaverage data (template data)\n",
    "overlap=np.zeros_like(thresholds)\n",
    "surf=os.path.join(pat_subjects_dir,'fsaverage_sym','surf','lh.inflated')\n",
    "surface=mt.load_mesh_geometry(surf)\n",
    "areas=os.path.join(pat_subjects_dir,'fsaverage_sym','surf','lh.area')\n",
    "area=mt.load_mesh_data(areas)\n",
    "neighbours=get_neighbours_from_tris(surface['faces'])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create empty variables\n",
    "# This variable will list whether there is an overlap between the cluster \n",
    "#and lesion label (i.e. correctly detected by classifer)\n",
    "any_overlap=np.zeros(len(thresholds))\n",
    "#Lists whether there is an overlap between the first (top) cluster \n",
    "#and lesion label (i.e. correctly detected by classifer)\n",
    "first_cluster_overlap=np.zeros(len(thresholds))\n",
    "#Lists number of predicted clusters\n",
    "number_of_clusters=np.zeros(len(thresholds))\n",
    "\n",
    "#Hemispheres\n",
    "hemis=['rh','lh']\n",
    "# For every patient\n",
    "for subject_id in patient_ids:\n",
    "    print(subject_id)\n",
    "    # For every threshold\n",
    "    for k, threshold in enumerate(thresholds):\n",
    "        clusters_stacked=[]\n",
    "        lesions_stacked=[]\n",
    "        predictions_stacked=[]\n",
    "        # for each hemisphere\n",
    "        for hemi in hemis:\n",
    "            # Change path and overlay to neural network output\n",
    "            test_nn_output=os.path.join(pat_subjects_dir,subject_id,'xhemi','classifier',hemi+'.'+subject_id+'.classifier_output.mgh')\n",
    "            \n",
    "            # Change path and file name to lesion label\n",
    "            mask=os.path.join(pat_subjects_dir,subject_id,'xhemi','surf',hemi+'.lesion_on_lh.mgh')\n",
    "\n",
    "            prediction=io.load_mgh(test_nn_output)\n",
    "            #if no lesion mask set all to zero\n",
    "            try:\n",
    "                lesion_mask=io.load_mgh(mask)\n",
    "            except IOError:\n",
    "                lesion_mask=np.zeros_like(prediction)\n",
    "            #clusters above threshold\n",
    "            clusters=cluster_above_threshold(prediction, np.array(neighbours), threshold=threshold)\n",
    "            #removes clusters below minimum area (e.g. 50mm)\n",
    "            clusters_filtered=filter_cluster(clusters,area,min_area_threshold)\n",
    "            #within a subject creates a variable containing all clusters across hemispheres\n",
    "            clusters_stacked.extend(clusters_filtered)\n",
    "            #creates variable with lesion across both hemispheres\n",
    "            lesions_stacked.extend(lesion_mask)\n",
    "            #variable of neural network output across both hemis\n",
    "            predictions_stacked.extend(prediction)\n",
    "            #turns previous variables (lists) into arrays\n",
    "        predictions_stacked=np.array(predictions_stacked)\n",
    "        lesions_stacked=np.array(lesions_stacked)\n",
    "        clusters_stacked=np.array(clusters_stacked)\n",
    "        #ranks clusters\n",
    "        clusters_ranked=rank_clusters(clusters_stacked,predictions_stacked,ranking=\"skew\")\n",
    "        #finds number of clusters\n",
    "        number_of_clusters[k]=np.max(clusters_ranked)\n",
    "        #finds if any clusters overlap with lesion label (ie, lesion detected by classifier)\n",
    "        any_overlap[k]+=check_any_overlap_with_lesion_mask(clusters_ranked,lesions_stacked)\n",
    "        # finds if first cluster overlaps with lesion label\n",
    "        first_cluster_overlap[k]+=check_cluster_overlap_with_lesion_mask(clusters_ranked,lesions_stacked,cluster_index=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Gets specificity at every threshold\n",
    "cont_any_lesions=np.zeros(len(thresholds))\n",
    "\n",
    "hemis=['rh','lh']\n",
    "# For each control participant\n",
    "for subject_id in control_ids:\n",
    "    print(subject_id)\n",
    "    for k, threshold in enumerate(thresholds):\n",
    "        clusters_stacked=[]\n",
    "        predictions_stacked=[]\n",
    "        for hemi in hemis:\n",
    "            #Change path and file name to control classifier output\n",
    "            test_nn_output=os.path.join(cont_subjects_dir,subject_id,'xhemi','classifier',hemi+'.'+subject_id+'.classifier_output.mgh')\n",
    "            prediction=io.load_mgh(test_nn_output)\n",
    "            #Clusters classifier output above the threshold\n",
    "            clusters=cluster_above_threshold(prediction, np.array(neighbours), threshold=threshold)\n",
    "            #Removes clusters below the threshold\n",
    "            clusters_filtered=filter_cluster(clusters,area,min_area_threshold)\n",
    "            # Creates list of clusters across both hemis\n",
    "            clusters_stacked.extend(clusters_filtered)\n",
    "            # Creates list of neural network output across both hemis\n",
    "            predictions_stacked.extend(prediction)\n",
    "            #Converts lists to arrays\n",
    "        predictions_stacked=np.array(predictions_stacked)\n",
    "        clusters_stacked=np.array(clusters_stacked)\n",
    "        #ranks clusters\n",
    "        clusters_ranked=rank_clusters(clusters_stacked,predictions_stacked,ranking=\"max\")\n",
    "        #Checks if each control has any detected areas\n",
    "        cont_any_lesions[k]+=(clusters_ranked>0).any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sensitivity is a variable containing the sensitivity of the classifier at all tested thresholds\n",
    "sensitivity=any_overlap/len(patient_ids)\n",
    "#Specificity is a variable containing the specificity of the classifier at all tested thresholds\n",
    "specificity=1-cont_any_lesions/len(control_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot sensitivity and specificity at different thresholds\n",
    "# x axis = thresholds\n",
    "# y axis = sensitivity + specificity\n",
    "%matplotlib inline\n",
    "plt.plot(thresholds,sensitivity+(specificity))\n",
    "\n",
    "# calculates optimal threshold - peak of the youden index\n",
    "optimal_threshold=get_optimal_threshold(sensitivity,specificity,thresholds,resolution=10)\n",
    "print(optimal_threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For the optimal theshold, calculate the number of clusters in each patient and whether lesion detected\n",
    "import nibabel as nb\n",
    "#Sets threshold as optimal threshold\n",
    "threshold = optimal_threshold\n",
    "hemis=['rh','lh']\n",
    "pat_number_of_clusters=np.zeros(len(patient_ids))\n",
    "pat_any_overlap=np.zeros(len(patient_ids))\n",
    "pat_first_cluster_overlap=np.zeros(len(patient_ids))\n",
    "\n",
    "# In every subject\n",
    "for k,subject_id in enumerate(patient_ids):\n",
    "    print(subject_id)\n",
    "    clusters_stacked=[]\n",
    "    lesions_stacked=[]\n",
    "    predictions_stacked=[]\n",
    "    # In each hemisphere\n",
    "    for h,hemi in enumerate(hemis):\n",
    "        #Path to classifier output in xhemi space(range of values 0 to 1 where values closer to 1 are more likely to be lesional)\n",
    "        test_nn_output=os.path.join(pat_subjects_dir,subject_id,'xhemi','classifier',hemi+'.'+subject_id+'.classifier_output.mgh')\n",
    "        #Path to lesion mask in xhemi space (this is the lesion mask created manually in native space and registered to xhemi)\n",
    "        mask=os.path.join(pat_subjects_dir,subject_id,'xhemi','surf',hemi+'.lesion_on_lh.mgh')\n",
    "        # This loads the classifier output\n",
    "        prediction=io.load_mgh(test_nn_output)\n",
    "        #This checks if there is a lesion mask for the hemisphere - if no lesion mask sets all values of variable prediction to zero\n",
    "        try:\n",
    "            lesion_mask=io.load_mgh(mask)\n",
    "        except IOError:\n",
    "             lesion_mask=np.zeros_like(prediction)\n",
    "        #This filters the classifier output above a threshold and then clusters the values above the threshold\n",
    "        clusters=cluster_above_threshold(prediction, np.array(neighbours), threshold=threshold)\n",
    "        #This filters out any clusters below the min_area_threshold (originally set to be 50mm2)\n",
    "        clusters_filtered=filter_cluster(clusters,area,min_area_threshold)\n",
    "        #This ranks the clusters. Clusters can be ranked according to mean, max, skew etc. \n",
    "        clusters_ordered=rank_clusters(clusters_filtered,prediction,ranking=\"skew\")\n",
    "        # This joins (groups) clusters that are within a number of vertices (n_steps) together\n",
    "        grouped=clump_clusters(clusters_ordered,np.array(neighbours),n_steps=5)\n",
    "        # This fills in holes in clusters e.g. if all neighbours are in the cluster, vertex will become part of cluster\n",
    "        defragged=defrag_surface(grouped, np.array(neighbours),steps=5)\n",
    "        defragged[defragged>0]+=h*1000\n",
    "\n",
    "        #This sets the final clusters into a variable called clusters_stacked which includes the data from both hemispheres\n",
    "        clusters_stacked.extend(defragged)\n",
    "        #This creates a variable called lesions_stacked that contains the lesion mask across both hemispheres\n",
    "        lesions_stacked.extend(lesion_mask)\n",
    "        #This creates a variable predictions_stacked that contains the original classifier output for both hemispheres\n",
    "        predictions_stacked.extend(prediction)\n",
    "    #This changes the variables predictions_stacked, lesions_stacked, clusters_stacked from a list into an array so that numpy commands can be used\n",
    "    predictions_stacked=np.array(predictions_stacked)\n",
    "    lesions_stacked=np.array(lesions_stacked)\n",
    "    clusters_stacked=np.array(clusters_stacked)\n",
    "    #This reranks the clusters across both hemispheres. Clusters can be ranked according to mean, max, skew, etc.\n",
    "    clusters_ranked=rank_clusters(clusters_stacked,predictions_stacked,ranking=\"skew\")\n",
    "    #This calculates the total number of clusters for the patient\n",
    "    pat_number_of_clusters[k]=np.max(clusters_ranked)\n",
    "    print(pat_number_of_clusters[k])\n",
    "    #This calculates whether any of the clusters overlap with the lesion label (sensitivity)\n",
    "    pat_any_overlap[k]=check_any_overlap_with_lesion_mask(clusters_ranked,lesions_stacked)\n",
    "    print(pat_any_overlap[k])\n",
    "    #This calculates whether the top cluster overlaps with the lesion label \n",
    "    pat_first_cluster_overlap[k]=return_cluster_overlap_with_lesion_mask(clusters_ranked,lesions_stacked)[0]\n",
    "    #This prints the cluster number that overlaps with the lesion label\n",
    "    print(pat_first_cluster_overlap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prints sensitivity, percentage of patients where the top cluster overlaps the lesion label\n",
    "# and prints average number of clusters per training patient at optimal threshold\n",
    "Sensitivity = np.mean(pat_any_overlap>0)\n",
    "print(Sensitivity)\n",
    "Percentage_first_cluster = np.mean(pat_first_cluster_overlap==1)\n",
    "print(Percentage_first_cluster)\n",
    "Mean_no_clusters =np.mean(pat_number_of_clusters)\n",
    "print(Mean_no_clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIND OUT SPECICITY IN CONTROLS\n",
    "#Create variable called control number of clusters\n",
    "cont_number_of_clusters=np.zeros(len(control_ids))\n",
    "\n",
    "# In every control\n",
    "for k,subject_id in enumerate(control_ids):\n",
    "    print(subject_id)\n",
    "    clusters_stacked=[]\n",
    "    predictions_stacked=[]\n",
    "    # In each hemisphere\n",
    "    for hemi in hemis:\n",
    "        #Path to classifier output in xhemi space(range of values 0 to 1 where values closer to 1 are more likely to be lesional)\n",
    "        test_nn_output=os.path.join(cont_subjects_dir,subject_id,'xhemi','classifier',hemi+'.'+subject_id+'.classifier_output.mgh')\n",
    "        # This loads the classifier output\n",
    "        prediction=io.load_mgh(test_nn_output)\n",
    "        #This filters the classifier output above a threshold and then clusters the values above the threshold\n",
    "        clusters=cluster_above_threshold(prediction, np.array(neighbours), threshold=threshold)\n",
    "        #This filters out any clusters below the min_area_threshold (originally set to be 50mm2)\n",
    "        clusters_filtered=filter_cluster(clusters,area,min_area_threshold)\n",
    "        #This ranks the clusters. Clusters can be ranked according to mean, max, skew etc. \n",
    "        clusters_ordered=rank_clusters(clusters_filtered,prediction,ranking=\"skew\")\n",
    "        # This joins (groups) clusters that are within a number of vertices (n_steps) together\n",
    "        grouped=clump_clusters(clusters_ordered,np.array(neighbours),n_steps=5)\n",
    "        # This fills in holes in clusters e.g. if all neighbours are in the cluster, vertex will become part of cluster\n",
    "        defragged=defrag_surface(grouped, np.array(neighbours),steps=5)\n",
    "        defragged[defragged>0]+=k*1000\n",
    "\n",
    "        #This sets the final clusters into a variable called clusters_stacked which includes the data from both hemispheres\n",
    "        clusters_stacked.extend(defragged)\n",
    "        #This creates a variable predictions_stacked that contains the original classifier output for both hemispheres\n",
    "        predictions_stacked.extend(prediction)\n",
    "    #This changes the variables predictions_stacked, lesions_stacked, clusters_stacked from a list into an array so that numpy commands can be used\n",
    "    predictions_stacked=np.array(predictions_stacked)\n",
    "    clusters_stacked=np.array(clusters_stacked)\n",
    "    #This reranks the clusters across both hemispheres. Clusters can be ranked according to mean, max, skew, etc.\n",
    "    clusters_ranked=rank_clusters(clusters_stacked,predictions_stacked,ranking=\"skew\")\n",
    "    #This calculates the total number of clusters for the patient\n",
    "    cont_number_of_clusters[k]=np.max(clusters_ranked)\n",
    "    print(cont_number_of_clusters[k])\n",
    "\n",
    "#Prints the number of clusters in each control participant\n",
    "print(cont_number_of_clusters)\n",
    "Mean_no_clusters_controls =np.mean(cont_number_of_clusters)\n",
    "#Prints the average number of clusters per control\n",
    "print(Mean_no_clusters_controls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nibabel as nb\n",
    "def process_predictions_to_ranked_clusters(seeg_subjects_dir,subject_id, neighbours, threshold, ranking=\"skew\",steps=5):\n",
    "    \"\"\"Load in classifier output (range of values between 0 and 1 where values closer to 1 more likely to be lesional)\n",
    "     Filter, cluster and rank the clusters across both hemispheres\n",
    "     Save out ranked clusters as\"\"\"\n",
    "    clusters_stacked=[]\n",
    "    predictions_stacked=[]\n",
    "    # In each hemisphere\n",
    "    hemis=['rh','lh']\n",
    "    for k, hemi in enumerate(hemis):\n",
    "        #Path to classifier output in xhemi space(range of values 0 to 1 where values closer to 1 are more likely to be lesional)\n",
    "        #test_nn_output=os.path.join(seeg_subjects_dir,subject_id,'xhemi','classifier',hemi+'.'+subject_id+'.classifier_output.mgh')\n",
    "        test_nn_output=os.path.join(seeg_subjects_dir,subject_id,'xhemi','classifier',hemi+'.'+subject_id+'.classifier_ouput.mgh')\n",
    "        # This loads the classifier output\n",
    "        prediction=io.load_mgh(test_nn_output)\n",
    "        #This filters the classifier output above a threshold and then clusters the values above the threshold\n",
    "        clusters=cluster_above_threshold(prediction, np.array(neighbours), threshold=threshold)\n",
    "        #This filters out any clusters below the min_area_threshold (originally set to be 50mm2)\n",
    "        clusters_filtered=filter_cluster(clusters,area,min_area_threshold)\n",
    "        #This ranks the clusters. Clusters can be ranked according to mean, max, skew etc. \n",
    "        clusters_ordered=rank_clusters(clusters_filtered,prediction,ranking=ranking)\n",
    "        print(clusters_ordered)\n",
    "        # This joins (groups) clusters that are within a number of vertices (n_steps) together\n",
    "        grouped=clump_clusters(clusters_ordered,np.array(neighbours),n_steps=steps)\n",
    "        # This fills in holes in clusters e.g. if all neighbours are in the cluster, vertex will become part of cluster\n",
    "\n",
    "        defragged=defrag_surface(grouped, np.array(neighbours),steps=steps)\n",
    "        defragged[defragged>0]+=k*1000\n",
    "\n",
    "        #This sets the final clusters into a variable called clusters_stacked which includes the data from both hemispheres\n",
    "        clusters_stacked.extend(defragged)\n",
    "        #This creates a variable predictions_stacked that contains the original classifier output for both hemispheres\n",
    "        predictions_stacked.extend(prediction)\n",
    "    #This changes the variables predictions_stacked, lesions_stacked, clusters_stacked from a list into an array so that numpy commands can be used\n",
    "    predictions_stacked=np.array(predictions_stacked)\n",
    "    clusters_stacked=np.array(clusters_stacked)\n",
    "    print(np.unique(clusters_stacked))\n",
    "    #This reranks the clusters across both hemispheres. Clusters can be ranked according to mean, max, skew, etc.\n",
    "    clusters_ranked=rank_clusters(clusters_stacked,predictions_stacked,ranking=ranking)\n",
    "    clusters_ranked_rh=clusters_ranked[:len(clusters_ranked)/2]\n",
    "    clusters_ranked_lh=clusters_ranked[len(clusters_ranked)/2:]\n",
    "    \n",
    "    \n",
    "    print(\"right\",np.unique(clusters_ranked_rh))\n",
    "    print(\"left\",np.unique(clusters_ranked_lh))\n",
    "    print(\"both\", np.unique(clusters_ranked))\n",
    "    demo=nb.load(test_nn_output)\n",
    "    # Saves an overlay for each hemi of the clusters after thresholding\n",
    "    io.save_mgh(os.path.join(seeg_subjects_dir,subject_id,'xhemi','classifier','lh.'+subject_id+'.classifier_output_clusters_after_thresholding.mgh'),clusters_ranked_lh,demo)\n",
    "    io.save_mgh(os.path.join(seeg_subjects_dir,subject_id,'xhemi','classifier','rh.'+subject_id+'.classifier_output_clusters_after_thresholding.mgh'),clusters_ranked_rh,demo)\n",
    "    return "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in SEEG data, threshold, cluster, plot electrodes and plot clusters and save as html and finally evaluate proximity of clusters to ictal contacts\n",
    "\n",
    "excel_path='/path/to/excel'\n",
    "## Load in your excel file containing seeg ictal, interictal, spread fiducials\n",
    "xls = pd.ExcelFile(os.path.join(excel_path,'sEEG_data.xlsx'))\n",
    "## Change path to location of sEEG data \n",
    "meld_path = '/path/to/seeg'\n",
    "## Change to location of freesurfer data for sEEG subjects\n",
    "freesurfer_dir = '/path/to/Freesurfer/'\n",
    "## List any patients you want to exclude\n",
    "exclude=np.array([['FCD_patient_1'],['FCD_patient_2']])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in interictal, ictal, all fiducials and affected hemisphere from excel spreadsheet\n",
    "Interictal = pd.read_excel( xls, 'Interictal' )\n",
    "Ictal = pd.read_excel( xls, 'Ictal')\n",
    "Fiducials = pd.read_excel( xls, 'Fiducials')\n",
    "Hemi = pd.read_excel( xls, 'Hemi')\n",
    "Included = pd.read_excel( xls, 'Included')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get subject_ids for seeg patients\n",
    "subject_ids =Included.loc[:,'ID'].values\n",
    "\n",
    "# Set cluster number to assess proximity of ictal fiducials from \n",
    "cluster_no=1\n",
    "\n",
    "print(subject_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#preload fsaverage data\n",
    "overlap=np.zeros_like(thresholds)\n",
    "surf=os.path.join(pat_subjects_dir,'fsaverage_sym','surf','lh.inflated')\n",
    "surface=mt.load_mesh_geometry(surf)\n",
    "areas=os.path.join(pat_subjects_dir,'fsaverage_sym','surf','lh.area')\n",
    "area=mt.load_mesh_data(areas)\n",
    "neighbours=get_neighbours_from_tris(surface['faces'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in classifier output (range of values between 0 and 1 where values closer to 1 more likely to be lesional)\n",
    "# Filter, cluster and rank the clusters across both hemispheres\n",
    "# Save out ranked clusters as \"classifier_output_clusters_after_thresholding.mgh\"\n",
    "\n",
    "for sub in subject_ids:\n",
    "    if sub in exclude:\n",
    "        continue\n",
    "    print(sub)\n",
    "    process_predictions_to_ranked_clusters(freesurfer_dir,sub, neighbours, optimal_threshold)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use cell below to plot \n",
    "label_colors=np.array(['#aaaaaa', '#FF5733','#FFC300'])\n",
    "\n",
    "for sub in subject_ids:\n",
    "    #if os.path.isdir(os.path.join(meld_path,sub)):\n",
    "    if os.path.isdir(os.path.join(meld_path,sub+'_recon')):\n",
    "        print(sub)\n",
    "        if sub in exclude:\n",
    "            continue\n",
    "        #Gets affected hemisphere of patient from excel    \n",
    "        hemi = np.array(Hemi['Hemi'][Hemi['ID'] == sub ])[0]\n",
    "        #Gets list of fiducials (electrodes) for patient from excel\n",
    "        fiducials = get_fiducials(sub)\n",
    "        #Prints list of electrodes implanted\n",
    "        print(fiducials)\n",
    "        # Gets list of ictal contacts\n",
    "        ictals = get_electrodes(Ictal, sub)\n",
    "        #Gets list of interictal contacts\n",
    "        interictals = get_electrodes(Interictal, sub)\n",
    "        colours=[]\n",
    "        electrode_coords=[]\n",
    "        for f in fiducials:\n",
    "            #Loads .fcsv file of electrode\n",
    "            electrode = np.genfromtxt(os.path.join(meld_path, sub, 'Slicer',f+'.fcsv'),delimiter=',', usecols = (0,1,2,3,4))\n",
    "            #check if only one electrode\n",
    "            if np.ndim(electrode)==1:\n",
    "                length_seeg = 1\n",
    "                #finds electrode contact-co-ordinates\n",
    "                electrode_coords=electrode[1:4]\n",
    "            else:\n",
    "                length_seeg=len(electrode[:,1:4])\n",
    "                if len(electrode_coords)==0:\n",
    "                    electrode_coords=electrode[:,1:4]\n",
    "                else:\n",
    "                    #Finds electrode contact co-ords\n",
    "                    electrode_coords=np.vstack((electrode_coords,electrode[:,1:4]))\n",
    "            #Codes ictal contacts red, interictals yellow and other contacts black\n",
    "            for k in range(length_seeg):\n",
    "                contact_name = f + str(k+1)\n",
    "                if contact_name in ictals:\n",
    "                    contact_value = 'red'\n",
    "                    print(contact_name)\n",
    "                elif contact_name in interictals:\n",
    "                    contact_value = 'yellow'\n",
    "                   # print(contact_name)\n",
    "                else :\n",
    "                    contact_value = 'black'\n",
    "                colours.append(contact_value)\n",
    "\n",
    "## The next part of the script creates an interactive plot and saves as an .html of the electrodes and the classifier output\n",
    "#It uses nilearn plotting.\n",
    "        connectome_info = _get_markers(electrode_coords, colours)\n",
    "        connectome_info[\"marker_size\"] = 10 #marker_size\n",
    "        plot_info = {\"connectome\": connectome_info}\n",
    "\n",
    "        hemis=['right', 'left']\n",
    "        hemifs=['rh','lh']\n",
    "        replacement='[\"left\", \"right\",\"labelleft\",\"labelright\"]'\n",
    "\n",
    "\n",
    "\n",
    "        for h,hemi in enumerate(hemis):\n",
    "        #plot gray translucent surface\n",
    "            pialname=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.pial')\n",
    "            whitename=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.white')\n",
    "\n",
    "\n",
    "            lesion_name=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.'+sub+'.classifier_output_clusters_after_thresholding.mgh')\n",
    "            lesion_surf=os.path.join(meld_path,sub+'_recon','surf',hemifs[h]+'.lesions.pial')\n",
    "            plot_info['pial_'+hemi]=mesh_to_plotly(pialname)    \n",
    "            neighbours=get_neighbours_from_tris(io_mesh.load_mesh_geometry(pialname)['faces'])\n",
    "            surf_map=tidy_integer_labels(np.array(io.import_mgh(lesion_name)),neighbours)\n",
    "\n",
    "            #set lesions above 1 to all same colour (2)\n",
    "            surf_map[surf_map>1]=2\n",
    "            plot_info['vertexcolor_'+hemi] = list(label_colors[surf_map])\n",
    "            plot_info['opacity_'+hemi] = 0.1\n",
    "\n",
    "            #plot opaque lesion labels\n",
    "            label=np.where(surf_map!=0)[0]  \n",
    "            label2object(lesion_surf,pialname,whitename,label)\n",
    "            # try because if no lesions on one side, no mesh is created and html fails\n",
    "            #so remove labelleft or labelright to plot\n",
    "            try:\n",
    "                plot_info['pial_label'+hemi] = mesh_to_plotly(lesion_surf)\n",
    "                plot_info['vertexcolor_label'+hemi] = list(np.hstack((label_colors[surf_map],label_colors[surf_map])))\n",
    "                plot_info['opacity_label'+hemi] = 0.4\n",
    "            except ValueError:\n",
    "                if hemi=='right':\n",
    "                    replacement='[\"left\", \"right\",\"labelleft\"]'\n",
    "                elif hemi=='left':\n",
    "                    replacement='[\"left\", \"right\",\"labelright\"]'\n",
    "\n",
    "        as_json = json.dumps(plot_info)\n",
    "        as_html = get_html_template('connectome_plot_template.html')\n",
    "        as_html= as_html.replace('[\"left\", \"right\"]',replacement)\n",
    "\n",
    "        as_html= as_html.replace('info[\"color\"] = \"#aaaaaa\";','info[\"vertexcolor\"] = connectomeInfo[\"vertexcolor_\" + hemisphere];')\n",
    "        as_html= as_html.replace('info[\"opacity\"] = getOpacity();','info[\"opacity\"] = connectomeInfo[\"opacity_\" + hemisphere];')\n",
    "\n",
    "        as_html=as_html.replace(\n",
    "            'INSERT_CONNECTOME_JSON_HERE', as_json)\n",
    "        as_html = add_js_lib(as_html, embed_js=True)\n",
    "\n",
    "        class ConnectomeView(HTMLDocument):\n",
    "            pass\n",
    "        ConnectomeView(as_html).save_as_html(os.path.join(meld_path,sub,sub+'.html'))\n",
    "        \n",
    "## The next part of the script plots a raincloud plot of distance of the ictal, interictal \n",
    "#and other contacts from a cluster. This is plotted for clusters 1 to 5.\n",
    "\n",
    "        for counter in range(5):\n",
    "            #not interested in number 0\n",
    "            cluster_no=counter+1\n",
    "            lesion_coords=[[0]]\n",
    "            for h,hemi in enumerate(hemis):\n",
    "                pialname=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.pial')\n",
    "                whitename=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.white')\n",
    "                #lesion_name=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.'+sub+'.classifier_output_clusters_after_thresholding.mgh')\n",
    "                lesion_name=os.path.join(meld_path,sub+'_recon','surf/'+hemifs[h]+'.'+sub+'.classifier_output_clusters_after_thresholding.mgh')\n",
    "                dum=lesion_coordinates(pialname,whitename,lesion_name,cluster_no)\n",
    "                if len(dum) >0:\n",
    "                    lesion_coords=dum\n",
    "            if lesion_coords[0][0]==0:\n",
    "                print(\"no cluster found for this cluster number, exiting loop\")\n",
    "                break\n",
    "            distances=np.zeros(np.shape(colours))\n",
    "            for k,coord in enumerate(electrode_coords):\n",
    "                closest_coord, distances[k]=closest_node(coord,lesion_coords)\n",
    "            classification=[]\n",
    "            poss_colours=['black','yellow','red']\n",
    "            for node in colours:\n",
    "                classification.append(poss_colours.index(node))\n",
    "            f, ax = plt.subplots(figsize=(12, 11))\n",
    "            colors=['black','orange','red']\n",
    "            plt.rcParams.update({'font.size': 30})\n",
    "            threshold=10\n",
    "            border=3\n",
    "            a=plt.plot([-2,3],[threshold,threshold],'r')\n",
    "            a=plt.fill_between([-2,3],[threshold-border,threshold-border],[threshold+border,threshold+border],facecolor='#ffdab9',alpha=0.2)\n",
    "\n",
    "            dy = \"Distance from predicted lesion (mm)\"; dx = \"Electrode classification\"; ort = \"v\"\n",
    "            ax=pt.half_violinplot(classification,distances, palette=sns.color_palette(colors), bw=.2,  linewidth=1,cut=0.,\\\n",
    "                               scale=\"area\", width=1, inner=None,orient=\"v\")\n",
    "            ax=sns.stripplot(classification,distances, palette=sns.color_palette(colors), edgecolor=\"white\",size=7,orient=\"v\",\\\n",
    "                             jitter=1,zorder=0)\n",
    "            ax=sns.boxplot(classification,distances, color=\"black\",orient=\"v\",width=.15,zorder=10,\\\n",
    "                          showcaps=True,boxprops={'facecolor':'none', \"zorder\":10},\\\n",
    "                           showfliers=True,whiskerprops={'linewidth':2, \"zorder\":10},saturation=1)\n",
    "\n",
    "            a=plt.ylabel(dy)\n",
    "            a=plt.xlabel(dx)\n",
    "\n",
    "            a=plt.xticks([0,1,2],['Nondescript','Interictal','Ictal'])\n",
    "            a=plt.savefig(os.path.join(meld_path,sub,sub+'_raincloud_'+str(cluster_no)+'.png'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}