Virtual_Lesion_Step1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Virtual Lesion\n",
    "This series of notebooks is used to determine the evidence for a pathway connecting two ROIs following those described in Pestilli et al.<sup>[1]</sup> and Leong et al.<sup>[2]</sup>.\n",
    "\n",
    "Required steps before running this notebook:\n",
    "\n",
    "<ol>\n",
    "    <li>Regions of interest (ROIs) were defined in MNI space and subsequently warped from MNI space to each individual subject’s native space</li>\n",
    "    <li>Constrained spherical deconvolution (CSD)<sup>[3]</sup> was used to estimate the fiber orientation distribution function (fODF). The fODFs were combined with tissue segmentations derived from the high resolution T1-weighted images to perform anatomically-constrained tractography<sup>[4]</sup></li>\n",
    "    <li>A whole-brain connectome was then generated by seeding all voxels in the white matter. These were filtered down using spherical-deconvolution informed filtering of tractograms (SIFT)<sup>[5]</sup> </li>\n",
    "</ol>\n",
    "\n",
    "Calculation steps done in this notebook:\n",
    "\n",
    "<ol>\n",
    "    <li>Generate fiber tracks seeding from ROI1 that pass through ROI2 and vise versa with MRtrix3</li>\n",
    "    <li>Remove the fibers from the previous step that do not start in one ROI and end in the other (The fibers that pass through the ROIs but do not terminate in them)</li>\n",
    "    <li>Combine all valid fiber tracks into one candiadate streamline set</li>\n",
    "    <li>Cluster the candiadate streamlines with dipys<sup>[6]</sup> QuickBundles algorithm<sup>[7]</sup> to remove outliers</li>\n",
    "</ol>\n",
    "\n",
    "<sup>[1]</sup> <i>PMID: 25194848</i> <br/>\n",
    "<sup>[2]</sup> <i>PMID: 26748088</i> <br/>\n",
    "<sup>[3]</sup> <i>PMID: 18583153</i> <br/>\n",
    "<sup>[4]</sup> <i>PMID: 22705374</i> <br/>\n",
    "<sup>[5]</sup> <i>PMID: 23238430</i> <br/>\n",
    "<sup>[6]</sup> <i>PMID: 24600385</i> <br/>\n",
    "<sup>[7]</sup> <i>PMID: 23248578</i> <br/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from nibabel import trackvis as tv\n",
    "from dipy.segment.clustering import QuickBundles\n",
    "import utilities\n",
    "from dipy.segment.metric import ResampleFeature\n",
    "from dipy.segment.metric import AveragePointwiseEuclideanMetric\n",
    "import itertools\n",
    "import MRTrix2TrackVis\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:interface:Reading header data...\n",
      "INFO:interface:...adding \"/data/hcp/data/100307/5TT.mif\" to header for key \"act\"\n",
      "INFO:interface:...adding \"1\" to header for key \"backtrack\"\n",
      "INFO:interface:...adding \"1\" to header for key \"crop_at_gmwmi\"\n",
      "INFO:interface:...adding \"3\" to header for key \"downsample_factor\"\n",
      "INFO:interface:...adding \"0.25\" to header for key \"fod_power\"\n",
      "INFO:interface:...adding \"0.1\" to header for key \"init_threshold\"\n",
      "INFO:interface:...adding \"8\" to header for key \"lmax\"\n",
      "INFO:interface:...adding \"45\" to header for key \"max_angle\"\n",
      "INFO:interface:...adding \"250\" to header for key \"max_dist\"\n",
      "INFO:interface:...adding \"500000000\" to header for key \"max_num_attempts\"\n",
      "INFO:interface:...adding \"2500\" to header for key \"max_num_tracks\"\n",
      "INFO:interface:...adding \"1\" to header for key \"max_seed_attempts\"\n",
      "INFO:interface:...adding \"1000\" to header for key \"max_trials\"\n",
      "INFO:interface:...adding \"iFOD2\" to header for key \"method\"\n",
      "INFO:interface:...adding \"2.5\" to header for key \"min_dist\"\n",
      "INFO:interface:...adding \"0.3.12-1060-gdbdc1293\" to header for key \"mrtrix_version\"\n",
      "INFO:interface:...adding \"0.625\" to header for key \"output_step_size\"\n",
      "INFO:interface:...adding \"0\" to header for key \"rk4\"\n",
      "INFO:interface:...adding \"4\" to header for key \"samples_per_step\"\n",
      "INFO:interface:...adding \"1\" to header for key \"sh_precomputed\"\n",
      "INFO:interface:...adding \"/data/hcp/data/100307/FOD.mif\" to header for key \"source\"\n",
      "INFO:interface:...adding \"0.625\" to header for key \"step_size\"\n",
      "INFO:interface:...adding \"0\" to header for key \"stop_on_all_include\"\n",
      "INFO:interface:...adding \"0.1\" to header for key \"threshold\"\n",
      "INFO:interface:...adding \"1458009798.6861047745\" to header for key \"timestamp\"\n",
      "INFO:interface:...adding \"0\" to header for key \"unidirectional\"\n",
      "INFO:interface:...adding \"Float32LE\" to header for key \"datatype\"\n",
      "INFO:interface:...adding \". 628\" to header for key \"file\"\n",
      "INFO:interface:...adding \"803\" to header for key \"count\"\n",
      "INFO:interface:...adding \"803\" to header for key \"total_count\"\n",
      "INFO:interface:Reached the end of the header!\n",
      "INFO:interface:Identifying the number of points per tract...\n",
      "INFO:interface:MRTrix Header:\n",
      "INFO:interface:{'mrtrix_version': '0.3.12-1060-gdbdc1293', 'backtrack': '1', 'max_dist': '250', 'file': '. 628', 'threshold': '0.1', 'output_step_size': '0.625', 'max_num_attempts': '500000000', 'downsample_factor': '3', 'source': '/data/hcp/data/100307/FOD.mif', 'step_size': '0.625', 'method': 'iFOD2', 'rk4': '0', 'min_dist': '2.5', 'init_threshold': '0.1', 'timestamp': '1458009798.6861047745', 'samples_per_step': '4', 'crop_at_gmwmi': '1', 'offset': 628, 'fod_power': '0.25', 'stop_on_all_include': '0', 'max_angle': '45', 'count': 803, 'max_seed_attempts': '1', 'max_num_tracks': '2500', 'sh_precomputed': '1', 'total_count': '803', 'max_trials': '1000', 'datatype': 'Float32LE', 'lmax': '8', 'unidirectional': '0', 'act': '/data/hcp/data/100307/5TT.mif'}\n",
      "INFO:interface:Applying transformation from scanner coordinates to /hcp/100307/T1w/Diffusion/data.nii.gz\n",
      "INFO:interface:Reading tracks...\n",
      "INFO:interface:100% : 803 tracks read\n",
      "INFO:interface:Saving Trackvis file as /data/hcp/data/100307/RVLPFC2FIRSTamyg_bigRight_combined.trk\n",
      "INFO:interface:TrackVis Header:\n",
      "INFO:interface:('TRACK', [145, 174, 145], [1.25, 1.25, 1.25], [0.0, 0.0, 0.0], 0, ['', '', '', '', '', '', '', '', '', ''], 0, ['', '', '', '', '', '', '', '', '', ''], [[-1.25, 0.0, 0.0, 90.0], [0.0, 1.25, 0.0, -126.0], [0.0, 0.0, 1.25, -72.0], [0.0, 0.0, 0.0, 1.0]], '', 'LAS', '', [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], '', '', '', '', '', '', '', 803, 2, 1000)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Process subject 100307\n",
      "    Clustered File does not exist for this subject, start calculation.\n",
      "    All neccessary files there, continue ...\n",
      "    Can access raw nifti data, start conversion and clustering.\n",
      "    Convert MRTrix streams to TrackVis\n",
      "    Cluster Steams\n",
      "    All done\n"
     ]
    }
   ],
   "source": [
    "# Library of Files\n",
    "path = '/hcp/'\n",
    "path_saveing = '/data/hcp/data/'\n",
    "\n",
    "subjects = os.listdir(path_saveing)\n",
    "\n",
    "subjects_sorted = sorted(subjects)\n",
    "subjects_sorted.remove('.nii.gz')\n",
    "\n",
    "for subject in subjects_sorted:\n",
    "    print 'Process subject ' + subject\n",
    "    \n",
    "    if os.path.isfile(os.path.join(path_saveing, subject, 'RVLPFC2FIRSTamyg_bigRight_clustered2.trk')) == False:\n",
    "        print \"    Clustered File does not exist for this subject, start calculation.\"\n",
    "    \n",
    "        if os.path.isfile(os.path.join(path_saveing, subject, 'FOD.mif')) == True and os.path.isfile(os.path.join(path_saveing, subject, 'ROI_RVLPFC_15mm_54_27_12.nii.gz')) == True and os.path.isfile(os.path.join(path_saveing, subject, 'ROI_FIRSTamyg_bigRight.nii.gz')) == True:\n",
    "            print \"    All neccessary files there, continue ...\"\n",
    "    \n",
    "            directory_output = os.path.join(path_saveing, subject)\n",
    "\n",
    "            if os.path.isfile(os.path.join(path_saveing, subject, 'RVLPFC2FIRSTamyg_bigRight_combined.tck')) == False:\n",
    "                print '    Fiber Tracks do not exist, start First Fiber Fracking'\n",
    "                cmd = \"tckgen \" + directory_output + \"/FOD.mif \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight1.tck -number 2500 -seed_image \" + directory_output + \"/ROI_RVLPFC_15mm_54_27_12.nii.gz  -include \" + directory_output + \"/ROI_FIRSTamyg_bigRight.nii.gz -force -maxnum 500000000 -act \" + directory_output + \"/5TT.mif -backtrack -crop_at_gmwmi -maxlength 250\"\n",
    "                os.system(cmd)\n",
    "\n",
    "                print '    Start Second Fiber Fracking'\n",
    "                cmd = \"tckgen \" + directory_output + \"/FOD.mif \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight2.tck -number 2500 -seed_image \" + directory_output + \"/ROI_FIRSTamyg_bigRight.nii.gz -include \" + directory_output + \"/ROI_RVLPFC_15mm_54_27_12.nii.gz -force -maxnum 500000000 -act \" + directory_output + \"/5TT.mif -backtrack -crop_at_gmwmi -maxlength 250\"\n",
    "                os.system(cmd)\n",
    "\n",
    "                print '    First step to remove too long fiber from the first streamlines'\n",
    "                cmd = \"tckedit \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight1.tck \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight1_cut.tck -include \" + directory_output + \"/ROI_FIRSTamyg_bigRight.nii.gz -test_ends_only -force\"\n",
    "                os.system(cmd)\n",
    "\n",
    "                print '    Second step to remove too long fiber from the first streamlines'\n",
    "                cmd = \"tckedit \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight1_cut.tck \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight1_cut_cut.tck -include \" + directory_output + \"/ROI_RVLPFC_15mm_54_27_12.nii.gz -test_ends_only -force\"\n",
    "                os.system(cmd)\n",
    "\n",
    "                print '    First step to remove too long fiber from the second streamlines'\n",
    "                cmd = \"tckedit \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight2.tck \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight2_cut.tck -include \" + directory_output + \"/ROI_RVLPFC_15mm_54_27_12.nii.gz -test_ends_only -force\"\n",
    "                os.system(cmd)\n",
    "\n",
    "                print '    Second step to remove too long fiber from the second streamlines'\n",
    "                cmd = \"tckedit \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight2_cut.tck \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight2_cut_cut.tck -include \" + directory_output + \"/ROI_FIRSTamyg_bigRight.nii.gz -test_ends_only -force\"\n",
    "                os.system(cmd)\n",
    "\n",
    "                print '    Combine resulting streamlines'\n",
    "                cmd = \"tckedit \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight1_cut_cut.tck \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight2_cut_cut.tck \" + directory_output + \"/RVLPFC2FIRSTamyg_bigRight_combined.tck  -force\"\n",
    "                os.system(cmd)\n",
    "                \n",
    "            else:\n",
    "                f_in_nifti = os.path.join(path, subject, 'T1w/Diffusion/data.nii.gz')\n",
    "                f_in_stream = os.path.join(directory_output, 'RVLPFC2FIRSTamyg_bigRight_combined.tck')\n",
    "                f_out_converted = os.path.join(directory_output, 'RVLPFC2FIRSTamyg_bigRight_combined.trk')\n",
    "                f_out_clustered = os.path.join(directory_output, 'RVLPFC2FIRSTamyg_bigRight_clustered.trk')\n",
    "                f_out_centroids = os.path.join(directory_output, 'RVLPFC2FIRSTamyg_bigRight_centroids.trk')\n",
    "\n",
    "                if os.path.isfile(f_in_nifti) == True:\n",
    "                    print \"    Can access raw nifti data, start conversion and clustering.\"\n",
    "\n",
    "                    print '    Convert MRTrix streams to TrackVis'\n",
    "                    try: \n",
    "                        MRTrix2TrackVis.convert_tck2trk(f_in_stream, f_in_nifti, f_out_converted)\n",
    "                    except:\n",
    "                        print 'Could not convert .tck to .trk'\n",
    "\n",
    "                    print '    Cluster Steams'\n",
    "                    try: \n",
    "                        streams, hdr = tv.read(f_out_converted)\n",
    "                        streamlines = [i[0] for i in streams]\n",
    "\n",
    "                        feature = ResampleFeature(nb_points=50)\n",
    "                        metric = AveragePointwiseEuclideanMetric(feature=feature)\n",
    "                        qb = QuickBundles(threshold=10., metric=metric)\n",
    "                        clusters = qb.cluster(streamlines)\n",
    "\n",
    "                        major_cluster = clusters > 60\n",
    "                        major_path = []\n",
    "                        for j in range(len(clusters)):\n",
    "                            if major_cluster[j] == True:\n",
    "                                major_path.append([streamlines[i] for i in clusters[j].indices])\n",
    "                        major_streams = list(itertools.chain(*major_path))\n",
    "\n",
    "                        strm = ((sl, None, None) for sl in major_streams)\n",
    "                        tv.write(f_out_clustered, strm,  hdr_mapping=hdr)\n",
    "                        \n",
    "                        print '    All done'\n",
    "                        \n",
    "                    except:\n",
    "                        print '    Could not Cluster streams'\n",
    "                else:\n",
    "                    print \"    Could not load raw diffusion data, skip conversion and clustering.\"\n",
    "        else:\n",
    "            print \"    Some input files are missing, skip this subject.\"\n",
    "    else:\n",
    "        print \"    Clustered File exists already for this subject, skip calculation.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}