Common functions shared between ProbTrack notebooks


In [2]:
import os
from string import Template 
from datetime import datetime

import numpy as np
import nibabel

def processed_dir(subject, fs='fs0'):
    return os.path.join('/{}/New_Server/RCV/MR_Processed/{}'.format(fs, subject))

def sc_dir(subject, fs='fs0'):
    return os.path.join(processed_dir(subject, fs), 'StructConn')

def hardi_dir(subject):
    return os.path.join(processed_dir(subject), 'HARDI')

def sc_script_dir(subject):
    return os.path.join(sc_dir(subject), 'scripts')

def bedpost_script_path(subject):
    return os.path.join(sc_script_dir(subject), 'bedpostx.sh')

def recon_dir(subject):
    return os.path.join(processed_dir(subject), 'Recon')

def scratch_dir(subject):
    return sc_dir(subject, fs='scratch/burnsss1')

def seed_file_path(subject):
    return os.path.join(scratch_dir(subject), 'seeds.txt')

def get_setup_data(subject):
    replace = {}
    replace['sub_subject'] = subject
    replace['sub_bdir'] = sc_dir(subject)
    replace['sub_recondir'] = recon_dir(subject)
    replace['sub_date'] = datetime.now()
    replace['sub_seedlist'] = 'seeds.txt'
    return replace

label_order = ['G_rectus', 'G_subcallosal', 'S_suborbital', 'S_orbital_med-olfact',
            'G_orbital', 'S_orbital-H_Shaped', 'G_and_S_transv_frontopol', 'G_and_S_cingul-Ant',
            'G_and_S_frontomargin', 'G_front_sup', 'S_front_sup', 'S_front_middle',
            'G_front_middle', 'S_orbital_lateral', 'S_front_inf', 'G_front_inf-Triangul',
            'G_front_inf-Orbital', 'Lat_Fis-ant-Horizont', 'Lat_Fis-ant-Vertical', 'S_circular_insula_ant',
            'S_circular_insula_sup', 'G_insular_short', 'G_Ins_lg_and_S_cent_ins', 'S_precentral-inf-part',
            'G_front_inf-Opercular', 'S_precentral-sup-part', 'G_precentral', 'S_central',
            'G_and_S_subcentral', 'Lat_Fis-post', 'S_circular_insula_inf', 'G_temp_sup-Plan_polar',
            'Pole_temporal', 'G_temp_sup-G_T_transv', 'S_temporal_transverse', 'G_temp_sup-Plan_tempo',
            'G_temp_sup-Lateral', 'S_temporal_sup', 'G_temporal_middle', 'S_temporal_inf',
            'G_temporal_inf', 'S_collat_transv_ant', 'S_collat_transv_post', 'G_oc-temp_med-Parahip',
            'S_oc-temp_lat', 'G_oc-temp_lat-fusifor', 'G_pariet_inf-Supramar', 'G_postcentral',
            'S_postcentral', 'G_and_S_paracentral', 'G_parietal_sup', 'S_intrapariet_and_P_trans',
            'S_interm_prim-Jensen', 'G_pariet_inf-Angular', 'S_occipital_ant', 'G_and_S_occipital_inf',
            'G_occipital_middle', 'S_oc_sup_and_transversal', 'G_occipital_sup', 'G_cuneus',
            'S_oc_middle_and_Lunatus', 'Pole_occipital', 'S_oc-temp_med_and_Lingual', 'G_oc-temp_med-Lingual',
            'S_calcarine', 'S_parieto_occipital', 'G_precuneus', 'S_subparietal',
            'G_cingul-Post-dorsal', 'G_cingul-Post-ventral', 'S_pericallosal', 'S_cingul-Marginalis',
            'G_and_S_cingul-Mid-Post', 'G_and_S_cingul-Mid-Ant']

variables = Template("""bdir=${sub_bdir}
recon_dir=${sub_recondir}
fs_mgz=$recon_dir/mri/orig.mgz
str_mgz=$recon_dir/mri/rawavg.mgz
aparcaseg_mgz=${sub_recondir}/mri/aparc+aseg.mgz
fs=anat/fs.nii.gz
str=anat/str.nii.gz
str_bet=anat/str_bet.nii.gz
aparcaseg=anat/aparc+aseg.nii.gz
ventricles=anat/ventricles.nii.gz
lh_wm=anat/wm.lh.nii.gz
rh_wm=anat/wm.rh.nii.gz
waypoints=waypoints.txt
fa=dtifit/dtifit_FA.nii.gz

subject=${sub_subject}

# transforms
fs2str=bedpostx/xfms/fs2str.mat
str2fs=bedpostx/xfms/str2fs.mat
fa2fs=bedpostx/xfms/fa2fs.mat
fs2fa=bedpostx/xfms/fs2fa.mat
fa2str=bedpostx/xfms/fa2str.mat
str2fa=bedpostx/xfms/str2fa.mat

# files for seeds
tmp_area=area.txt
seed_list=${sub_seedlist}""")

def collapse_probtrack_results(waytotal_file, matrix_nii):
    with open(waytotal_file) as f:
        waytotal = int(f.read())
    data = nibabel.load(matrix_nii).get_data()    
    collapsed = data.sum(axis=0) / waytotal * 100.
    if waytotal == 0:
        collapsed = np.zeros(collapsed.shape)
    return collapsed

def single_process(subject_sc):
    template = os.path.join(subject_sc, 'results/{roi}.nii.gz.probtrackx2/matrix_seeds_to_all_targets.nii.gz')
    seeds_file = os.path.join(subject_sc, 'seeds.txt')
    f = open(seeds_file)
    with f:    
        processed_seed_list = [s.replace('.nii.gz','').replace('label/', '') 
                               for s in f.read().split('\n')
                               if s]
    N = len(processed_seed_list)
    conn = np.zeros((N, N))
    for idx, roi in enumerate(processed_seed_list):
        result = template.format(roi=roi)
        seed_directory = os.path.dirname(result)
        try:
            waytotal_file = os.path.join(seed_directory, 'waytotal')
            matrix_nii = os.path.join(seed_directory, 'matrix_seeds_to_all_targets.nii.gz')
            collapsed = collapse_probtrack_results(waytotal_file, matrix_nii)
            conn[idx, :] = collapsed
        except IOError:
            pass
    return conn, processed_seed_list, N

subject_list = ['061_206924',
                '063_207046',
                '064_207264',
                '067_207215',
                '072_207335',
#                 '073_207351', # bad reg
#                 '081_207486', # bad reg
#                 '093_207987', # bad reg
#                 '098_208247', # bad reg
#                 '106_208335', # bad reg
                '130_208994',
                '131_209154',
                '140_209143',
                '141_209157',
                '144_209407',
                '146_209355',
                '147_209378',
                '148_209625',
                '162_210032',
                '170_210044',
                '172_209736',
                '188_210443',
                '191_210512',
                '196_210780',
                '197_210808',
#                 '199_210894',
                '203_211015',
                '208_211122',
                '216_211291',
                '228_211662',]

subject_list.sort()

In [2]:
assert processed_dir('foo') == '/fs0/New_Server/RCV/MR_Processed/foo'
assert sc_dir('foo') == '/fs0/New_Server/RCV/MR_Processed/foo/StructConn'
assert hardi_dir('foo') == '/fs0/New_Server/RCV/MR_Processed/foo/HARDI'
assert sc_script_dir('foo') == '/fs0/New_Server/RCV/MR_Processed/foo/StructConn/scripts'
assert bedpost_script_path('foo') == '/fs0/New_Server/RCV/MR_Processed/foo/StructConn/scripts/bedpostx.sh'
assert recon_dir('foo') == '/fs0/New_Server/RCV/MR_Processed/foo/Recon'
assert len(label_order) == 74