Merge pull request #280 from berkeley-stat159/master

merge from master to final

code/utils/scripts/multi_beta_script.py

@@ -0,0 +1,96 @@
+
+"""
+Purpose:
+-----------------------------------------------------------------------------------
+We generate beta values for each single voxels for each subject and save them to 
+files for multi-comparison test.
+-----------------------------------------------------------------------------------
+
+"""
+
+
+import sys, os
+##sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
+sys.path.append(os.path.join(os.path.dirname('__file__'), "../functions/"))
+import numpy as np
+from glm import *
+#from convolution_normal_script import X_matrix
+#from convolution_high_res_script import X_matrix_high_res
+import nibabel as nib
+import matplotlib.pyplot as plt
+from smoothing import *
+from t_test import *
+
+# Create the necessary directories if they do not exist
+dirs = ['../../../txt_output/multi_beta']
+for d in dirs:
+    if not os.path.exists(d):
+            os.makedirs(d)
+
+# Locate the different paths
+project_path = '../../../'
+# TODO: change it to relevant path
+conv_path = project_path + 'txt_output/conv_normal/'
+conv_high_res_path = project_path + 'txt_output/conv_high_res/'
+
+# select your own subject
+subject_list = [str(i) for i in range(1,17)]
+
+conv_list = [str(i) for i in range(1,5)]
+
+txt_paths = [('ds005_sub' + s.zfill(3) + '_t1r1' +'_cond'+ c.zfill(3),\
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv001_canonical.txt', \
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv002_canonical.txt', \
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv003_canonical.txt', \
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv004_canonical.txt', \
+              '../../../data/ds005_2/sub' + s.zfill(3) + '/model/model001/task001_run002' \
+              #for me (Min) it's ds005_2... the filtered data set
+              + '.feat/filtered_func_data_mni.nii.gz',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv_001_high_res.txt',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv_002_high_res.txt',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv_003_high_res.txt',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r2' +'_conv_004_high_res.txt') \
+                for s in subject_list \
+                for c in conv_list]
+
+for txt_path in txt_paths:
+# get 4_d image data
+    name = txt_path[0]
+
+    img = nib.load(txt_path[5])
+    data = img.get_data()
+
+    p = 7
+    # p is the number of columns in our design matrix
+    # it is the number of convolved column plus 1 (a column of 1's)
+
+    X_matrix1 = np.loadtxt(txt_path[1])
+    X_matrix2 = np.loadtxt(txt_path[2])
+    X_matrix3 = np.loadtxt(txt_path[3])
+    X_matrix4 = np.loadtxt(txt_path[4])
+    X_matrix = np.ones((len(X_matrix1),p))
+    X_matrix[...,1] = X_matrix1
+    X_matrix[...,2] = X_matrix2
+    X_matrix[...,3] = X_matrix3
+    X_matrix[...,4] = X_matrix4
+    linear_drift = np.linspace(-1, 1, 240)
+    quadratic_drift = linear_drift ** 2
+    quadratic_drift -= np.mean(quadratic_drift)
+    X_matrix[...,5] = linear_drift
+    X_matrix[...,6] = quadratic_drift
+
+    # smooth the data
+    # use high resolution matrix and re-run the regression
+    data_smooth = smoothing(data,1,range(data.shape[-1]))
+    beta_3d_smooth, t, df, p = t_stat(data_smooth,X_matrix)
+    beta_3d_smooth_task = beta_3d_smooth[...,1]
+    beta_3d_smooth_gain = beta_3d_smooth[...,2]
+    beta_3d_smooth_loss = beta_3d_smooth[...,3]
+    beta_3d_smooth_dist = beta_3d_smooth[...,4]
+
+    location_of_txt= dirs[0]
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_task.txt",beta_3d_smooth_task.ravel())
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_gain.txt",beta_3d_smooth_gain.ravel())
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_loss.txt",beta_3d_smooth_loss.ravel())
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_dist.txt",beta_3d_smooth_dist.ravel())
+

code/utils/scripts/multi_comparison_script.py

@@ -0,0 +1,245 @@
+
+"""
+Purpose:
+-----------------------------------------------------------------------------------
+We seek the activated voxel positionsi through multi-comparison of beta values across
+subjects
+
+Step
+-----------------------------------------------------------------------------------
+1. calculate the mean of each single beta values across subject and plot them
+2. calculate the variance of each single beta values across subject and plot them
+3. calculate the t-stat of each single beta values across subject and plot them
+4. calculate the p-value of each single betav values across subject and plot them
+"""
+
+
+import sys, os
+##sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
+sys.path.append(os.path.join(os.path.dirname('__file__'), "../functions/"))
+sys.path.append(os.path.join(os.path.dirname('__file__'), "./"))
+import numpy as np
+from glm import *
+import nibabel as nib
+from matplotlib import colors
+import matplotlib.pyplot as plt
+#from scipy.stats import sem
+from smoothing import *
+from plot_mosaic import * 
+from scipy.stats import t as t_dist
+#from visualization import *
+
+
+nice_cmap_values = np.loadtxt('actc.txt')
+nice_cmap = colors.ListedColormap(nice_cmap_values, 'actc')
+dirs = ['../../../txt_output/multi_beta','../../../fig/multi_beta']
+plt.rcParams['image.cmap'] = 'gray'
+plt.rcParams['image.interpolation'] = 'nearest'
+project_path='../../../'
+
+for d in dirs:
+    if not os.path.exists(d):
+        os.makedirs(d)
+
+
+#need for plotting mask
+#TODO: change the location path
+template = nib.load(project_path+\
+           'data/mni_icbm152_t1_tal_nlin_asym_09c_2mm.nii')
+template_data = template.get_data()
+img = nib.load(project_path+'data/ds005_2/sub001/model/model001/task001_run002.feat/' + \
+ 'masked_filtered_func_data_mni.nii.gz')
+
+task = dict()
+gain = dict()
+loss = dict()
+
+#load all of them
+for x in range(1,17):
+	task[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_task.txt')
+
+for x in range(1,17):
+	gain[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_gain.txt')
+
+for x in range(1,17):
+	loss[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_loss.txt')
+
+# for x in range(1,17):
+# 	dist[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_dist.txt')
+
+
+##################################### MEAN plot #########################################
+
+#calculate mean and plot (let's try for task)
+task_sum = task[1]
+for x in range(2,17):
+	task_sum +=task[x]
+
+task_mean = task_sum/16
+task_mean_reshape = task_mean.reshape(91,109,91)
+
+
+gain_sum = gain[1]
+for x in range(2,17):
+	gain_sum +=gain[x]
+
+gain_mean = gain_sum/16
+gain_mean_reshape = gain_mean.reshape(91,109,91)
+
+
+loss_sum = gain[1]
+for x in range(2,17):
+	loss_sum +=loss[x]
+
+loss_mean = loss_sum/16
+loss_mean_reshape = loss_mean.reshape(91,109,91)
+
+
+plt.title('In brain activated voxels - \nmean across 16 subjects on TASK condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(task_mean_reshape, transpose=False), cmap='seismic', alpha=1, vmin=task_mean_reshape.min(), vmax= task_mean_reshape.max())
+plt.colorbar()
+plt.savefig(dirs[1]+'/mean_task.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nmean across 16 subjects on GAIN condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(gain_mean_reshape, transpose=False), cmap='seismic', alpha=1, vmin=gain_mean_reshape.min(), vmax= gain_mean_reshape.max())
+plt.colorbar()
+plt.savefig(dirs[1]+'/mean_gain.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nmean across 16 subjects on LOSS condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(loss_mean_reshape, transpose=False), cmap='seismic', alpha=1, vmin=loss_mean_reshape.min(), vmax= loss_mean_reshape.max())
+plt.colorbar()
+plt.savefig(dirs[1]+'/mean_loss.png')
+plt.clf()
+
+
+##################################### SD plot #########################################
+
+#calculate variance and plot
+stdlst = []
+for x in range(1,17):
+	stdlst.append(task[x])
+
+stdarray = np.array(stdlst)
+task_std = stdarray.std(axis=0)
+#task_std.shape -> (902629,0)
+task_std_reshape = task_std.reshape(91,109,91)
+
+stdlst = []
+for x in range(1,17):
+	stdlst.append(gain[x])
+
+stdarray = np.array(stdlst)
+gain_std = stdarray.std(axis=0)
+#task_std.shape -> (902629,0)
+gain_std_reshape = gain_std.reshape(91,109,91)
+
+stdlst = []
+for x in range(1,17):
+	stdlst.append(loss[x])
+
+stdarray = np.array(stdlst)
+loss_std = stdarray.std(axis=0)
+#task_std.shape -> (902629,0)
+loss_std_reshape = loss_std.reshape(91,109,91)
+
+plt.title('In brain activated voxels - \nStandard Deviation across 16 subjects on TASK condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(task_std_reshape, transpose=False), cmap='seismic', alpha=1, vmin=task_std_reshape.min(), vmax= task_std_reshape.max())
+plt.colorbar()
+plt.savefig(dirs[1]+'/std_task.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nStandard Deviation across 16 subjects on GAIN condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(gain_std_reshape, transpose=False), cmap='seismic', alpha=1, vmin=gain_std_reshape.min(), vmax= gain_std_reshape.max())
+plt.colorbar()
+plt.savefig(dirs[1]+'/std_gain.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nStandard Deviation across 16 subjects on LOSS condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(loss_std_reshape, transpose=False), cmap='seismic', alpha=1, vmin=loss_std_reshape.min(), vmax= loss_std_reshape.max())
+plt.colorbar()
+plt.savefig(dirs[1]+'/std_loss.png')
+plt.clf()
+
+
+##################################### stat plot #########################################
+
+
+#calculate t-stat and plot
+task_tstat = task_mean/(task_std/np.sqrt(15))
+task_tstat = np.nan_to_num(task_tstat)
+task_tstat_reshape = task_tstat.reshape(91,109,91)
+
+gain_tstat = gain_mean/(gain_std/np.sqrt(15))
+gain_tstat = np.nan_to_num(gain_tstat)
+gain_tstat_reshape = gain_tstat.reshape(91,109,91)
+
+loss_tstat = loss_mean/(loss_std/np.sqrt(15))
+loss_tstat = np.nan_to_num(loss_tstat)
+loss_tstat_reshape = loss_tstat.reshape(91,109,91)
+
+
+plt.title('In brain activated voxels - \nT-statistics across 16 subjects on TASK condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(task_tstat_reshape, transpose=False), cmap='seismic', alpha=1)
+plt.colorbar()
+plt.savefig(dirs[1]+'/tstat_task.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nT-statistics across 16 subjects on GAIN condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(gain_tstat_reshape, transpose=False), cmap='seismic', alpha=1)
+plt.colorbar()
+plt.savefig(dirs[1]+'/tstat_gain.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nT-statistics across 16 subjects on LOSS condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(loss_tstat_reshape, transpose=False), cmap='seismic', alpha=1)
+plt.colorbar()
+plt.savefig(dirs[1]+'/tstat_loss.png')
+plt.clf()
+
+
+
+##################################### P-value plot #########################################
+
+#calculate p-value and plot
+task_pval = t_dist.cdf(abs(task_tstat), 15)
+task_pval_reshape = task_pval.reshape(91,109,91)
+task_pval_reshape = 1-task_pval_reshape
+
+gain_pval = t_dist.cdf(abs(gain_tstat), 15)
+gain_pval_reshape = gain_pval.reshape(91,109,91)
+gain_pval_reshape = 1-gain_pval_reshape
+
+loss_pval = t_dist.cdf(abs(loss_tstat), 15)
+loss_pval_reshape = loss_pval.reshape(91,109,91)
+loss_pval_reshape = 1-loss_pval_reshape
+p_value_thres = 0.05/(91*109*91)
+
+plt.title('In brain activated voxels - \nP-value across 16 subjects on TASK condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(task_pval_reshape, transpose=False), cmap='Greys', alpha=1)
+plt.colorbar()
+plt.savefig(dirs[1]+'/pval_task.png')
+plt.clf()
+
+plt.title('In brain activated voxels - \nP-value across 16 subjects on GAIN condition', fontsize=12)
+plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)
+plt.imshow(plot_mosaic(gain_pval_reshape, transpose=False), cmap='Greys', alpha=1)
+plt.colorbar()
+plt.savefig(dirs[1]+'/pval_gain.png')
+plt.clf()
+
+
+
+
+