Skip to content

Calculate weighted average values for core skeleton and peripheral skeleton (outside of ROIs) in ENIGMA analysis #55

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CaitlinLloyd wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Calculate weighted average values for core skeleton and peripheral skeleton (outside of ROIs) in ENIGMA analysis #55

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3ce46d8
Delete roi_analysis.py
CaitlinLloyd May 3, 2023
221817d
Delete tbss_all
CaitlinLloyd May 3, 2023
0a8f911
Add files via upload
CaitlinLloyd May 3, 2023
219b47d
Update roi_analysis.py
CaitlinLloyd May 11, 2023
4a561d4
Update roi_analysis.py
CaitlinLloyd May 21, 2023
385a684
Update roi_analysis.py
CaitlinLloyd May 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 57 additions & 11 deletions lib/roi_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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)
Expand All @@ -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=[]
Expand All @@ -104,15 +114,16 @@ 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)


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())

Expand All @@ -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()
Expand All @@ -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)
Expand All @@ -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)
Expand Down
14 changes: 13 additions & 1 deletion lib/tbss_all
100755 → 100644
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand All @@ -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, '
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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)



Expand Down