#!/usr/bin/env python2 # -*- coding: utf-8 -*- ''' author: Christian Olaf Häusler created on Wednesday September 15 2021 start ipython2.7 via 'python2 -m IPython' ''' from __future__ import print_function from collections import OrderedDict from glob import glob from mvpa2.datasets.mri import fmri_dataset from nilearn.input_data import NiftiMasker, MultiNiftiMasker import argparse import numpy as np import os import pandas as pd import matplotlib.pyplot as plt import nibabel as nib import re import subprocess # from mvpa2.datasets.mri import fmri_dataset from scipy.stats import gaussian_kde # constants GRP_PPA_PTTRN = 'sub-??/masks/in_bold3Tp2/grp_PPA_bin.nii.gz' AO_FOV_MASK = 'sub-??/masks/in_bold3Tp2/audio_fov.nii.gz' GM_MASK = 'sub-??/masks/in_bold3Tp2/gm_bin_dil_fov.nii.gz' # binary mask(s) of individual visual localizer (in subject space) PPA_MASK_PATTERN = 'inputs/studyforrest-data-visualrois/'\ 'sub-*/rois/?PPA_?_mask.nii.gz' # individual 2nd level results (primary cope in subject space) PREDICTED_PREFIX = 'predicted-VIS-PPA_from_' PREDICTIONS = [ 'anatomy.nii.gz', 'srm-ao-av-vis_feat10_0-451.nii.gz', 'srm-ao-av-vis_feat10_0-892.nii.gz', 'srm-ao-av-vis_feat10_0-1330.nii.gz', 'srm-ao-av-vis_feat10_0-1818.nii.gz', 'srm-ao-av-vis_feat10_0-2280.nii.gz', 'srm-ao-av-vis_feat10_0-2719.nii.gz', 'srm-ao-av-vis_feat10_0-3261.nii.gz', 'srm-ao-av-vis_feat10_0-3524.nii.gz', 'srm-ao-av-vis_feat10_3524-3975.nii.gz', 'srm-ao-av-vis_feat10_3524-4416.nii.gz', 'srm-ao-av-vis_feat10_3524-4854.nii.gz', 'srm-ao-av-vis_feat10_3524-5342.nii.gz', 'srm-ao-av-vis_feat10_3524-5804.nii.gz', 'srm-ao-av-vis_feat10_3524-6243.nii.gz', 'srm-ao-av-vis_feat10_3524-6785.nii.gz', 'srm-ao-av-vis_feat10_3524-7123.nii.gz', 'srm-ao-av-vis_feat10_0-7123.nii.gz' ] VIS_ZMAP_PATTERN = 'inputs/studyforrest-data-visualrois/'\ 'sub-??/2ndlvl.gfeat/cope*.feat/stats/zstat1.nii.gz' # contrast used by Sengupta et al. (2016) to create the PPA mask VIS_VPN_COPES = OrderedDict({ # dicts are ordered from Python 3.7 'sub-01': 'cope8', 'sub-02': 'cope3', 'sub-03': 'cope3', 'sub-04': 'cope3', 'sub-05': 'cope3', 'sub-06': 'cope3', 'sub-09': 'cope3', 'sub-14': 'cope3', 'sub-15': 'cope3', 'sub-16': 'cope3', 'sub-17': 'cope3', 'sub-18': 'cope8', 'sub-19': 'cope3', 'sub-20': 'cope3' }) def parse_arguments(): ''' ''' parser = argparse.ArgumentParser( description='Creates mosaic of individual Bland-Altman-Plots' ) parser.add_argument('-i', default='./', help='input directory') parser.add_argument('-o', default='test', help='output directory') args = parser.parse_args() inDir = args.i outDir = args.o return inDir, outDir def find_files(pattern): ''' ''' found_files = glob(pattern) found_files = sort_nicely(found_files) return found_files def sort_nicely(l): '''Sorts a given list in the way that humans expect ''' convert = lambda text: int(text) if text.isdigit() else text alphanum_key = lambda key: [convert(c) for c in re.split('([0-9]+)', key)] l.sort(key=alphanum_key) return l def compute_means(data1, data2, log='n'): ''' ''' if len(data1) != len(data2): raise ValueError('data1 does not have the same length as data2.') data1 = np.asarray(data1) data2 = np.asarray(data2) if log == 'n': means = np.mean([data1, data2], axis=0) elif log == 'y': # what ever computation pass return means def compute_diffs(data1, data2, log='n'): ''' ''' if len(data1) != len(data2): raise ValueError('data1 does not have the same length as data2.') data1 = np.asarray(data1) data2 = np.asarray(data2) if log == 'n': diffs = data1 - data2 # Difference between data1 and data2 elif log == 'y': # what ever computation pass return diffs def process_df(subj, zmaps_df, out_path, prediction): ''' ''' # get the contrasts' names from the column names zmap1 = zmaps_df.iloc[:, 1] zmap2 = zmaps_df.iloc[:, 0] # mask all voxels not contained in the PPA group mask ppa_grp_masked1 = zmap1.as_matrix()[zmaps_df['ppa_grp'].as_matrix() > 0] ppa_grp_masked2 = zmap2.as_matrix()[zmaps_df['ppa_grp'].as_matrix() > 0] # mask all voxels not contained in the individual PPA mask ppa_ind_masked1 = zmap1.as_matrix()[zmaps_df['ppa_ind'].as_matrix() > 0] ppa_ind_masked2 = zmap2.as_matrix()[zmaps_df['ppa_ind'].as_matrix() > 0] # ao_ind_masked1 = zmap1.as_matrix()[zmaps_df['ao_ind'].as_matrix() > 0] # ao_ind_masked2 = zmap2.as_matrix()[zmaps_df['ao_ind'].as_matrix() > 0] datasets = [ [zmap1, zmap2], [ppa_grp_masked1, ppa_grp_masked2], [ppa_ind_masked1, ppa_ind_masked2], # [ao_ind_masked1, ao_ind_masked2] ] means_list = [compute_means(data1, data2) for data1, data2 in datasets] diffs_list = [compute_diffs(data1, data2) for data1, data2 in datasets] # Set up the axes with gridspec fig = plt.figure(figsize=(6, 6)) grid = plt.GridSpec(6, 6, hspace=0.0, wspace=0.0) # add three subplots ax_scatter = fig.add_subplot(grid[1:, :-1]) ax_xhist = fig.add_subplot(grid[0:1, 0:-1], yticklabels=[], sharex=ax_scatter) ax_yhist = fig.add_subplot(grid[1:, -1], xticklabels=[], sharey=ax_scatter) ax_scatter.text(5.1, 5.8, subj, fontsize=16, fontweight='bold') # plot voxel within occipitotemporal cortex plot_blandaltman(ax_scatter, means_list[0], diffs_list[0], alpha=0.6, c='darkgrey', s=2) # plot voxels within PPA group overlap plot_blandaltman(ax_scatter, means_list[1], diffs_list[1], alpha=1, c='royalblue', s=2) # plot voxels within individual PPA ROI plot_blandaltman(ax_scatter, means_list[2], diffs_list[2], alpha=1, c='r', s=2) # # plot voxels within (thresholded) individual AO zmap # plot_blandaltman(ax_scatter, # means_list[3], # diffs_list[3], # alpha=0.5, # c='y', # s=2) plot_histogram(ax_xhist, ax_yhist, means_list[0], diffs_list[0], alpha=1, color='darkgrey') plot_histogram(ax_xhist, ax_yhist, means_list[1], diffs_list[1], alpha=1, color='royalblue') plot_histogram(ax_xhist, ax_yhist, means_list[2], diffs_list[2], alpha=1, color='r') # try: # plot_histogram(ax_xhist, ax_yhist, # means_list[3], diffs_list[3], # alpha=1, # color='y') # # except ValueError: # print(subj, 'has no significant cluster in primary AO contrast') # save it suffix = 'from_' + prediction.split('.nii')[0] out_file = ('%s_bland-altman_%s.png' % (subj, suffix)) out_fpath = os.path.join(out_path, out_file) if not os.path.exists(out_path): os.makedirs(out_path) plt.savefig(out_fpath, bbox_inches='tight', dpi=80) plt.close() def plot_blandaltman(ax, means, diffs, *args, **kwargs): ''' ''' if len(means) != len(diffs): raise ValueError('means do not have the same length as diffs.') # annotation # variable subj is still a global here if subj in ['sub-01', 'sub-04', 'sub-09', 'sub-16', 'sub-19']: ax.set_ylabel('Difference between 2 measures', fontsize=16) if subj in ['sub-19', 'sub-20']: ax.set_xlabel('Average of 2 measures', fontsize=16) # draw the scattergram ax.scatter(means, diffs, *args, **kwargs) # set the size of the tick labels ax.xaxis.set_tick_params(labelsize=16) ax.yaxis.set_tick_params(labelsize=16) # limit the range of data shown ax.set_xlim(-5, 5) ax.set_ylim(-5, 5) # draw horizontal and vertical line at 0 ax.axhline(0, color='k', linewidth=0.5, linestyle='--') ax.axvline(0, color='k', linewidth=0.5, linestyle='--') return None def plot_histogram(x_hist, y_hist, means, diffs, *args, **kwargs): ''' ''' # basic preparation of the axes x_hist.xaxis.set_tick_params(bottom=True, labelbottom=False) y_hist.yaxis.set_tick_params(bottom=True, labelbottom=False) # show vertical/horizontal zero line x_hist.axvline(0, color='k', linewidth=0.5, linestyle='--') y_hist.axhline(0, color='k', linewidth=0.5, linestyle='--') # plot histogram -> take KDE plot, s. below # x_hist.hist(means, 50, normed=True, histtype='bar', # orientation='vertical', **kwargs) # y_hist.hist(diffs, 50, normed=True, histtype='bar', # orientation='horizontal', **kwargs) # x_hist.set_yscale('log') # y_hist.set_xscale('log') # plot KDE xvalues = np.arange(-5.1, 5.1, 0.1) kde_means = gaussian_kde(means) kde_diffs = gaussian_kde(diffs) # KDE subplot on the top x_hist.plot(xvalues, kde_means(xvalues), **kwargs) # x_hist.fill_between(xvalues, kde_means(xvalues), 0, **kwargs) x_hist.set_ylim(0.015, 0.8) x_hist.set_yticks([0.2, 0.4, 0.6]) # , 1]) x_hist.set_yticklabels(['.2', '.4', '.6']) # , '1']) x_hist.yaxis.set_tick_params(labelsize=16) # KDE subplot on the right y_hist.plot(kde_diffs(xvalues), xvalues, **kwargs) # y_hist.fill_between(xvalues, kde_diffs(xvalues), 0, **kwargs) y_hist.set_xlim(0.015, .8) y_hist.set_xticks([0.2, 0.4, 0.6]) # , 1]) y_hist.set_xticklabels(['.2', '.4', '.6']) # , '1']) y_hist.xaxis.set_tick_params(labelsize=16) return None def create_mosaic(in_pattern, dims, out_fpath, dpi): ''' http://www.imagemagick.org/Usage/montage/ ''' # dimensions in columns*rows subprocess.call( ['montage', '-density', str(dpi), in_pattern, '-geometry', '+1+1', '-tile', dims, out_fpath]) def load_mask(subj): ''' ''' # open the mask of cortices (at the moment it justs the union of # individual PPAs) grp_ppa_fpath = GRP_PPA_PTTRN.replace('sub-??', subj) # print('PPA GRP mask:\t', grp_ppa_fpath) # subject-specific field of view in audio-description ao_fov_mask = AO_FOV_MASK.replace('sub-??', subj) # print('ind. FoV mask:\t', ao_fov_mask) # load the masks grp_ppa_img = nib.load(grp_ppa_fpath) ao_fov_img = nib.load(ao_fov_mask) # (dilated) gray matter mask; see constat at script's top #gm_fpath = GM_MASK.replace('sub-??', subj) #gm_img = nib.load(gm_fpath) #final_mask_data = grp_ppa_img.get_fdata() * gm_img.get_fdata() final_mask_data = grp_ppa_img.get_fdata() * ao_fov_img.get_fdata() final_mask_img = nib.Nifti1Image(final_mask_data, grp_ppa_img.affine, header=grp_ppa_img.header) return final_mask_img if __name__ == "__main__": # get command line argument inDir, outDir = parse_arguments() vis_fpathes = find_files(VIS_ZMAP_PATTERN) # get pathes & filenames of all available zmaps subjs_pattern = os.path.join(inDir, 'sub-??') subjs_pathes = find_files(subjs_pattern) subjs = [re.search(r'sub-..', string).group() for string in subjs_pathes] subjs = sorted(list(set(subjs))) for prediction in PREDICTIONS[:]: # loop over subjects for subj in subjs: print('\nProcessing', subj) # initialize a dataframe zmaps_df = pd.DataFrame() # load predicted z-map of current subject pred_fpath = os.path.join(inDir, subj, PREDICTED_PREFIX + prediction) pred_fpath = pred_fpath.replace('sub-??', subj) pred_data = fmri_dataset(pred_fpath).samples.T # put the array into the dataframe zmaps_df['predicted'] = np.ndarray.flatten(pred_data) # load the visual localizer's z-map # filter for current subject fpathes = [x for x in vis_fpathes if subj in x] # filter for the subject's correct cope vis_fpath = [x for x in fpathes if VIS_VPN_COPES[subj] in x][0] vis_data = fmri_dataset(vis_fpath).samples.T # put the array into the dataframe zmaps_df['vis'] = np.ndarray.flatten(vis_data) # load the GRP PPA mask (in subject space) # (GRP PPA will be masked with individual FoV) final_mask_img = load_mask(subj) final_mask_data = final_mask_img.get_fdata() # put the array into the dataframe zmaps_df['ppa_grp'] = np.ndarray.flatten(final_mask_data) # load the subject-specific PPA mask (Sengupta et al., 2016) ppa_fpathes = find_files(PPA_MASK_PATTERN.replace('sub-*', subj)) ppaIndData = fmri_dataset(ppa_fpathes).samples.sum(axis=0) # masked it with individual FoV (just to be sure) ao_fov_fpath = find_files(AO_FOV_MASK.replace('sub-??', subj)) ao_fov_data = fmri_dataset(ao_fov_fpath).samples ppaIndMasked = ppaIndData * ao_fov_data zmaps_df['ppa_ind'] = np.ndarray.flatten(ppaIndMasked) # let the function do the calculation and the plotting process_df(subj, zmaps_df, outDir, prediction) # create the mosaic suffix = 'from_' + prediction.split('.nii')[0] infile_pattern = os.path.join(outDir, ('sub-??_bland-altman_%s.png' % suffix)) out_fpath = os.path.join(outDir, 'subjs_bland-altman_%s.png' % suffix) create_mosaic(infile_pattern, '3x5', out_fpath, 80) # columns x rows