diff --git a/lib/roi_analysis.py b/lib/roi_analysis.py index aa090be..0429748 100644 --- a/lib/roi_analysis.py +++ b/lib/roi_analysis.py @@ -48,7 +48,7 @@ def average_labels(labels): return np.unique(commonLabels) -def subject_stat(imgPath, c, modality, label2name, commonLabels, labelMap, roiDir, avgFlag): +def subject_stat(imgPath, c, modality, label2name, commonLabels, labelMap, roiDir, avgFlag,skeletonMaskPeri): print('Creating ROI based statistics for', imgPath) img= load(imgPath).get_fdata() @@ -71,7 +71,16 @@ def subject_stat(imgPath, c, modality, label2name, commonLabels, labelMap, roiDi df.loc[i+1]= [label2name[intLabel]]+ [num2str(x) for x in [_img_roi.mean(), _img_roi.size]] else: df.loc[i + 1] = [label2name[intLabel]] + ['0','0'] - + + roi= load(skeletonMaskPeri).get_fdata() + _roi = np.logical_and(_imgNonzero, roi) + _img_roi= img[_roi] + end=len(df) + if _img_roi.size: + df.loc[end]=["Peri"]+ [num2str(x) for x in [_img_roi.mean(), _img_roi.size]] + else: + df.loc[end] = ["Peri"] + ['0','0'] + df.set_index('Tract').to_csv(stat_file) # FIXME: save unsorted df to match with that of ENIGMA? # df.sort_values(by='Tract').set_index('Tract').to_csv(stat_file) @@ -82,8 +91,9 @@ def subject_stat(imgPath, c, modality, label2name, commonLabels, labelMap, roiDi df_avg = pd.DataFrame(columns=['Tract', 'Average', 'nVoxels']) df_avg.loc[0] = df.loc[0].copy() + df_avg.loc[1] = df.loc[end].copy() - row= 1 + row= 2 for common in commonLabels: dm=[] num=[] @@ -104,7 +114,8 @@ def subject_stat(imgPath, c, modality, label2name, commonLabels, labelMap, roiDi str(int(np.sum(num)))] row = row + 1 break - + + # FIXME: save unsorted df_avg so Tract, Right-Tract, and Left-Tract are together? df_avg.sort_values(by='Tract').set_index('Tract').to_csv(avg_stat_file) print('Made ', avg_stat_file) @@ -112,7 +123,7 @@ def subject_stat(imgPath, c, modality, label2name, commonLabels, labelMap, roiDi def roi_analysis(imgs, cases, args, roiDir, N_CPU): - intLabels = load(args.labelMap).get_fdata() + intLabels = load(args.labelMap).get_data() label2name = parse_labels(np.unique(intLabels)[1:], args.lut) commonLabels= average_labels(label2name.values()) @@ -121,7 +132,7 @@ def roi_analysis(imgs, cases, args, roiDir, N_CPU): # subject_stat(imgPath, c, args.modality, label2name, commonLabels, intLabels, roiDir, args.avg) pool.apply_async(func= subject_stat, args= (imgPath, c, args.modality, label2name, commonLabels, - intLabels, roiDir, args.avg), error_callback= RAISE) + intLabels, roiDir, args.avg, args.skeletonMaskPeri), error_callback= RAISE) pool.close() pool.join() @@ -132,12 +143,27 @@ def roi_analysis(imgs, cases, args, roiDir, N_CPU): # avg_stat_file= pjoin(roiDir, f'{c}_{modality}_roi_avg.csv') # read one stat_file, obtain headers df= pd.read_csv(pjoin(roiDir, f'{cases[0]}_{args.modality}_roi.csv')) - df_comb= pd.DataFrame(columns= np.append('Cases', df['Tract'].values)) - + df_comb= pd.DataFrame(columns= np.append(['Cases','Weighted_avg','Core_weighted_avg'], df['Tract'].values)) + print(df_comb) + for i, c in enumerate(cases): df= pd.read_csv(pjoin(roiDir, f'{c}_{args.modality}_roi.csv')) + core=df + core=core.drop(core[core['Tract'] =='Peri'].index, inplace = False) + core=core.drop(core[core['Tract'] =='AverageFA'].index, inplace = False) # num2str() text formatting is for precision control - df_comb.loc[i]= np.append(c, np.array([num2str(x) for x in df['Average'].values])) + no_avg=df + no_avg=no_avg.drop(no_avg[no_avg['Tract'] =='AverageFA'].index, inplace = False) + total_vox=sum(no_avg['nVoxels'].values) + no_avg['weight']=no_avg['Average'].values*(no_avg['nVoxels'].values/total_vox) + weighted_avg=sum(no_avg['weight'].values) + + + total_core=sum(core['nVoxels'].values) + core['weight']=core['Average'].values*(core['nVoxels'].values/total_core) + core_weighted_avg=sum(core['weight'].values) + print(np.append([c,weighted_avg,core_weighted_avg], np.array([num2str(x) for x in df['Average'].values]))) + df_comb.loc[i]= np.append([c,weighted_avg,core_weighted_avg], np.array([num2str(x) for x in df['Average'].values])) combined_stat= pjoin(args.statsDir, f'{args.modality}_combined_roi.csv') df_comb.sort_index(axis=1).set_index('Cases').to_csv(combined_stat) @@ -146,12 +172,32 @@ def roi_analysis(imgs, cases, args, roiDir, N_CPU): if args.avg: # read one avg_stat_file, obtain headers df_avg= pd.read_csv(pjoin(roiDir, f'{cases[0]}_{args.modality}_roi_avg.csv')) - df_avg_comb= pd.DataFrame(columns= np.append('Cases', df_avg['Tract'].values)) + df_avg_comb= pd.DataFrame(columns= np.append(['Cases','Weighted_avg','Core_weighted_avg'], df_avg['Tract'].values)) + + for i, c in enumerate(cases): + df = pd.read_csv(pjoin(roiDir, f'{c}_{args.modality}_roi_avg.csv')) + if args.avg: + # read one avg_stat_file, obtain headers + df_avg= pd.read_csv(pjoin(roiDir, f'{cases[0]}_{args.modality}_roi_avg.csv')) + df_avg_comb= pd.DataFrame(columns= np.append(['Cases','Weighted_avg','Core_weighted_avg'], df_avg['Tract'].values)) for i, c in enumerate(cases): df = pd.read_csv(pjoin(roiDir, f'{c}_{args.modality}_roi_avg.csv')) + core=df + core=core.drop(core[core['Tract'] =='Peri'].index, inplace = False) + core=core.drop(core[core['Tract'] =='AverageFA'].index, inplace = False) + # num2str() text formatting is for precision control + no_avg=df + no_avg=no_avg.drop(no_avg[no_avg['Tract'] =='AverageFA'].index, inplace = False) + total_vox=sum(no_avg['nVoxels'].values) + no_avg['weight']=no_avg['Average'].values*(no_avg['nVoxels'].values/total_vox) + weighted_avg=sum(no_avg['weight'].values) + total_core=sum(core['nVoxels'].values) + core['weight']=core['Average'].values*(core['nVoxels'].values/total_core) + core_weighted_avg=sum(core['weight'].values) + # num2str() text formatting is for precision control - df_avg_comb.loc[i] = np.append(c, np.array([num2str(x) for x in df['Average'].values])) + df_avg_comb.loc[i] = np.append([c,weighted_avg,core_weighted_avg], np.array([num2str(x) for x in df['Average'].values])) combined_avg_stat= pjoin(args.statsDir, f'{args.modality}_combined_roi_avg.csv') df_avg_comb.sort_index(axis=1).set_index('Cases').to_csv(combined_avg_stat) diff --git a/lib/tbss_all b/lib/tbss_all old mode 100755 new mode 100644 index 02dd47e..45ed9f6 --- a/lib/tbss_all +++ b/lib/tbss_all @@ -90,7 +90,7 @@ as you would provide for FA TBSS.""") help='create all of template, templateMask, skeleton, skeletonMask, and skeletonMaskDst') parser.add_argument('--enigma', action='store_true', - help='use ENGIMA provided template, templateMask, skeleton, skeletonMask, and skeletonMaskDst, ' + help='use ENGIMA provided template, templateMask, skeleton, skeletonMask, perimask and skeletonMaskDst, ' 'do JHU white matter atlas based ROI analysis using ENIGMA look up table') parser.add_argument('--fmrib', action='store_true', @@ -113,6 +113,9 @@ as you would provide for FA TBSS.""") parser.add_argument('--skeletonMaskDst', type=str, help='skeleton mask distance map') + + parser.add_argument('--skeletonMaskPeri', type=str, + help='skeleton mask peripheral') parser.add_argument('-s', '--space', type=str, help='you may register your template (including ANTs) to another standard space i.e MNI, ' @@ -199,6 +202,7 @@ as you would provide for FA TBSS.""") args.skeleton = None args.skeletonMask = None args.skeletonMaskDst = None + args.skeletonMaskPeri = None # obtain --caselist, --space, --labelMap, --lut, --avg, --SEARCH_RULE_MASK, --SKEL_THRESH from log # read the latest log @@ -280,6 +284,8 @@ as you would provide for FA TBSS.""") args.skeletonMask = abspath(args.skeletonMask) if args.skeletonMaskDst: args.skeletonMaskDst = abspath(args.skeletonMaskDst) + if args.skeletonMaskPeri: + args.skeletonMaskPeri = abspath(args.skeletonMaskPeri) if args.space: args.space = abspath(args.space) if args.labelMap: @@ -320,6 +326,7 @@ as you would provide for FA TBSS.""") args.templateMask = pjoin(enigmaDir, 'ENIGMA_DTI_FA_mask.nii.gz') args.skeleton = pjoin(enigmaDir, 'ENIGMA_DTI_FA_skeleton.nii.gz') args.skeletonMask = pjoin(enigmaDir, 'ENIGMA_DTI_FA_skeleton_mask.nii.gz') + args.skeletonMaskPeri = pjoin(enigmaDir, 'ENIGMA_DTI_FA_skeleton_mask_periphery.nii.gz') args.skeletonMaskDst = pjoin(enigmaDir, 'ENIGMA_DTI_FA_skeleton_mask_dst.nii.gz') args.lut = pjoin(enigmaDir, 'ENIGMA_look_up_table.txt') @@ -371,6 +378,11 @@ as you would provide for FA TBSS.""") fileBaseName= basename(args.skeletonMaskDst) if fileBaseName not in files: check_call(f'ln -s {args.skeletonMaskDst}', shell= True) + + if args.skeletonMaskPeri: + fileBaseName= basename(args.skeletonMaskPeri) + if fileBaseName not in files: + check_call(f'ln -s {args.skeletonMaskPeri}', shell= True)