kjamison / nemo

Network Modification Tool
MIT License
14 stars 3 forks source link

nemo_logo

NeMo 2.1 - Network Modification Tool

Predict brain network disruption from a lesion mask. Original concept described in Kuceyeski 2013.

NEW! Cloud interface for this tool can be used here: https://kuceyeski-wcm-web.s3.amazonaws.com/nemo/index.html

DOI

Contents

  1. Workflow overview
  2. Inputs
  3. Outputs
    1. Code examples for parsing outputs
  4. Website usage
  5. Requirements
  6. Details of tractography database
  7. Parcellations

Workflow overview

The general workflow for this tool consists of a database generation stage and a lesion disconnectivity stage:

  1. Tractography database generation: databaseprep/
    1. Compute whole-brain tractogram streamlines for 420 unrelated healthy subjects from the Human Connectome Project (See hcp_subjects_unrelated420_scfc.txt for list)
      • Both probabilistic (iFOD2+ACT) and deterministic (SD_STREAM) are available
    2. Nonlinearly warp streamlines into a common reference space (MNI152 v6, eg: $FSLDIR/data/standard/MNI152_T1_1mm.nii.gz, or website/atlases/MNI152_T1_1mm_brain.nii.gz)
    3. Additional processing to facillitate efficient computation: run_full_database_prep.sh
  2. Lesion to disconnectivity mapping: nemo_lesion_to_chaco.py
    1. Given a lesion mask, identify all streamlines it intersects, and identify the endpoints of those streamlines to compute brain regions for which we expect a reduction of connectivity
    2. Compute ChaCo (Change in Connectivity) score, which is the ratio of (disrupted streamlines)/(total streamlines) for each voxel or ROI (chacovol), or voxel/ROI pair (chacoconn). 0=no disconnection. 1=complete disconnection

We have created a user-friendly web interface to run this tool in the cloud (AWS):

Inputs

Outputs

Code examples for parsing outputs

Load chacovol_mean.pkl file and save as text/tsv/csv

import pickle
import numpy as np
data = pickle.load(open("mylesion_chacovol_mean.pkl","rb"))
#data contains a single ROIxROI matrix
np.savetxt("mylesion_chacovol_mean.txt",data,delimiter="\t")

#or save as .mat
from scipy.io import savemat
savemat("mylesion_chacovol_mean.mat",{"data":data})
# (load in matlab with M=load("mylesion_chacovol_mean.mat"); data=M.data;)

Load chacoconn_mean.pkl file and save as text/tsv/csv

#NOTE: chacoconn files are *SPARSE* format and must be converted to *DENSE* before saving to text using data.toarray()
import pickle
import numpy as np
data = pickle.load(open("mylesion_chacoconn_mean.pkl","rb"))
np.savetxt("mylesion_chacoconn_mean.txt",data.toarray(),delimiter="\t")

#or save dense/full matrix to matlab .mat:
from scipy.io import savemat
savemat("mylesion_chacoconn_mean.mat",{"data":data.toarray()})
# (load in matlab with M=load("mylesion_chacoconn_mean.mat"); data=M.data;)

#sparse chacoconn outputs *CAN* be saved as sparse format in .mat files *ONLY IF* you convert them to np.double:
from scipy.io import savemat
savemat("mylesion_chacoconn_mean_sparse.mat",{"data":data.astype(np.double)})
# (load in matlab with M=load("mylesion_chacoconn_mean_sparse.mat"); data=M.data; or data=full(M.data); ... )

Load chacovol_*_allref.pkl or chacovol_*_allref_denom.pkl files and save as matlab .mat

#These files contain SPARSE matrices, so we need to either convert them to dense:
import pickle
import numpy as np
from scipy.io import savemat
data = pickle.load(open("mylesion_chacovol_fs86subj_allref.pkl","rb"))
savemat("mylesion_chacovol_fs86subj_allref.mat",{"allref":data.toarray()})
# (load in matlab with M=load("mylesion_chacovol_fs86subj_allref.mat"); allref=M.allref;)

#or save as sparse after converting to np.double:
savemat("mylesion_chacovol_fs86subj_allref_sparse.mat",{"allref":data.astype(np.double)})
# (load in matlab with M=load("mylesion_chacovol_fs86subj_allref.mat"); allref=M.allref; or allref=full(M.allref); ...)

Load chacoconn_*_allref.pkl or chacoconn_*_allref_denom.pkl files and save as matlab .mat

#These files contain LISTS of SPARSE matrices, so we need to either
#a) convert them all to dense and stack them into 3D:
import pickle
import numpy as np
from scipy.io import savemat
data = pickle.load(open("mylesion_chacoconn_fs86subj_allref.pkl","rb"))
allref=np.stack([x.toarray() for x in data], axis=2) #creates an 86x86x420 3D matrix
savemat("mylesion_chacoconn_fs86subj_allref.mat",{"allref":allref})
# (load in matlab with M=load("mylesion_chacoconn_fs86subj_allref.mat"); allref=M.allref;)

#or b) save as a cell array of of sparse matrices to save space, but you MUST convert to np.double:
allref=[x.astype(np.double) for x in data]
savemat("mylesion_chacoconn_fs86subj_allref_sparse.mat",{"allref":allref})
# (load in matlab with M=load("mylesion_chacoconn_fs86subj_allref.mat"); allref=M.allref; allref_subj1=M.allref{1}; allref_subj1_full=full(M.allref{1}); ...)

Load chacoconn_resXmm_mean.pkl file for VOXEL-PAIRWISE output

#chacoconn_mean is a VOXELSxVOXELS pairwise matrix, where chacoconn_mean[i,:] represents 
# a VOLUME of disconnectivity to from voxel i to all other voxels
import pickle
import numpy as np
import nibabel as nib
from nilearn import plotting

#read the corresponding chacovol_resXmm_mean.nii.gz output which defines the output dimensions
Vref = nib.load("mylesion_chacovol_res5mm_mean.nii.gz")

data = pickle.load(open("mylesion_chacoconn_res5mm_mean.pkl","rb"))

#for this demonstration, find the voxel with the largest total disconnectivity, convert that row
#of the chacoconn output to a volume, and visualize that volume
data_sum=np.sum(data,axis=0)+np.sum(data,axis=1).T #data is upper triangular so we need to sum across both dimensions
voxel_index=np.argmax(data_sum)
newdata=data[voxel_index,:]+data[:,voxel_index].T
Vnew=nib.Nifti1Image(np.reshape(newdata.toarray(),Vref.shape),affine=Vref.affine, header=Vref.header)

#convert voxel index to coordinates we can include the position in the output filename
voxel_ijk=np.unravel_index(voxel_index,Vref.shape) #convert to (i,j,k) coords
voxel_xyz=(Vref.affine @ np.array(voxel_ijk+(1,)).T)[:3] #can also convert to x,y,z mm using Vref affine
voxel_ijk_string="%d_%d_%d" % (voxel_ijk)

#save glassbrain image of this volume using nilearn plotting tools
plotting.plot_glass_brain(Vnew,output_file="mylesion_chacoconn_res5mm_voxel_%s.png" % (voxel_ijk_string),colorbar=True)

#save volume as nii.gz so you can view it in nifti-viewing software (eg: fsleyes)
nib.save(Vnew,"mylesion_chacoconn_res5mm_voxel_%s.nii.gz" % (voxel_ijk_string))

#just for fun, save a second volume with JUST that seed voxel highlighted so we can overlay it later
cursordata=np.zeros(data.shape[0])
cursordata[voxel_index]=1
Vcursor=nib.Nifti1Image(np.reshape(cursordata,Vref.shape),affine=Vref.affine, header=Vref.header)
nib.save(Vcursor,"mylesion_chacoconn_res5mm_voxel_%s_cursor.nii.gz" % (voxel_ijk_string))

Website usage

Requirements

Details of tractography database

Parcellations

FS86 FS111 FS191 CocoMMP438 CocoMMPSUIT439 CocoYeo143 CocoYeo243 CocoYeo443 CocoLaus157 CocoLaus262 CocoLaus491 AAL AAL3 CC200 CC400 Shen268 Yeo7 Yeo17