Merge pull request #68 from LiamFengLin/liam_development

finished testing

code/network_analysis.py

@@ -9,7 +9,7 @@
 from gaussian_filter import spatial_smooth
 from general_utils import prepare_standard_img, prepare_mask, prepare_standard_data, form_cond_filepath
 from os.path import join
-from connectivity_utils import c_between, c_within
+from connectivity_utils import c_between, c_within, permute
 import numpy as np
 import os
 import math
@@ -161,64 +161,42 @@ def preprocessing_pipeline(subject_num, task_num, standard_source_prefix, cond_f
 
   return b_vols, img, in_brain_mask
 
-def subject_z_values(img, data, dist_from_center, dic, in_brain_mask):
+def subject_c_values(img, data, dist_from_center, dic, in_brain_mask):
 	mm_to_vox = npl.inv(img.affine)
 
 	roi_extractor = roi_extraction.SphereExtractor(in_brain_mask, dist_from_center)
 
 	expanded_dic = expand_dic(dic, mm_to_vox, roi_extractor)
 
-	mean_z_values = c_within(data, expanded_dic)
-	mean_z_values.update(c_between(data, expanded_dic))
-	return mean_z_values
+	mean_c_values = c_within(data, expanded_dic)
+	mean_c_values.update(c_between(data, expanded_dic))
+	return mean_c_values
 
-def group_z_values(standard_group_source_prefix, cond_filepath_prefix, dist_from_center, dic, group_info):
+def group_c_values(standard_group_source_prefix, cond_filepath_prefix, dist_from_center, dic, group_info):
 	task_nums = ("001", "002", "003")
 
 	# store layout
   #   level 1: task (0-back, 1-back, 2-back)
   #   level 2: group name (CON, SCZ)
   #   level 3: network name
   #   level 4: a list of ROI-ROI correlations
-	z_values_store = {"001":{"con":{}, "scz":{}},
+	c_values_store = {"001":{"con":{}, "scz":{}},
 										"002":{"con":{}, "scz":{}},
 										"003":{"con":{}, "scz":{}}}
 
 	for group, subject_nums in group_info.items():
 		for sn in subject_nums:
 			for tn in task_nums:
 				data, img, in_brain_mask = preprocessing_pipeline(sn, tn, standard_group_source_prefix, cond_filepath_prefix)
-				mean_z_values_per_net_pair = subject_z_values(img, data, dist_from_center, dic, in_brain_mask)
+				mean_c_values_per_net_pair = subject_c_values(img, data, dist_from_center, dic, in_brain_mask)
 
-				for network_pair_name, z_value in mean_z_values_per_net_pair.items():
+				for network_pair_name, c_value in mean_c_values_per_net_pair.items():
 					group_name = "con" if group in ("fmri_con", "fmri_con_sib") else "scz"
-					if network_pair_name not in z_values_store[tn][group_name]:
-						z_values_store[tn][group_name][network_pair_name] = [z_value]
+					if network_pair_name not in c_values_store[tn][group_name]:
+						c_values_store[tn][group_name][network_pair_name] = [c_value]
 					else:
-						z_values_store[tn][group_name][network_pair_name].append(z_value)
-	return z_values_store
-
-def permute (r1,r2):
-  """
-  This function performs the permuation test to two lists of r-values (r1:scz; r2:con).
-  Ho: mu_r1 = mu_r2
-  H1: muri < mu_r2
-  input:
-  1.r1 and r2 are two arrays containing the r-values
-  output:
-  1. one sided p-values
-  """
-  n1 = len(r1)
-  n2 = len(r2)
-  t_obs = np.mean(r1)-np.mean(r2)
-  pool = r1+r2
-  diff = []
-  for i in range(0,1000):
-    sample = random.sample(pool,n1)
-    diff.append(np.mean(sample)-(sum(pool)-sum(sample))/n2)
-  p_value = sum(list(i <= t_obs for i in diff))/len(diff)
-  return p_value
-
+						c_values_store[tn][group_name][network_pair_name].append(c_value)
+	return c_values_store
 
 if __name__ == "__main__":
 
@@ -236,28 +214,28 @@ def permute (r1,r2):
             "fmri_scz":("007", "009", "017", "031"),
             "fmri_scz_sib":("006", "008", "018", "024")}
 
-  z_values_store = group_z_values(standard_group_source_prefix, cond_filepath_prefix, dist_from_center, dic, small_group_info)
+  c_values_store = group_c_values(standard_group_source_prefix, cond_filepath_prefix, dist_from_center, dic, small_group_info)
 
-  generate_connectivity_results(z_values_store, output_filename)
+  generate_connectivity_results(c_values_store, output_filename)
 
   # change target r-values into list format
-  con_dmn_cer = np.ravel(z_values_store["003"]["con"]["Default-Cerebellar"]).tolist()
-  scz_dmn_cer = np.ravel(z_values_store["003"]["scz"]["Default-Cerebellar"]).tolist()
+  con_dmn_cer = np.ravel(c_values_store["003"]["con"]["Default-Cerebellar"]).tolist()
+  scz_dmn_cer = np.ravel(c_values_store["003"]["scz"]["Default-Cerebellar"]).tolist()
 
-  con_cer_co = np.ravel(z_values_store["003"]["con"]["Cerebellar-Cingulo-Opercular"]).tolist()
-  scz_cer_co = np.ravel(z_values_store["003"]["scz"]["Cerebellar-Cingulo-Opercular"]).tolist()
+  con_cer_co = np.ravel(c_values_store["003"]["con"]["Cerebellar-Cingulo-Opercular"]).tolist()
+  scz_cer_co = np.ravel(c_values_store["003"]["scz"]["Cerebellar-Cingulo-Opercular"]).tolist()
 
-  con_dmn_co = np.ravel(z_values_store["003"]["con"]["Default-Cingulo-Opercular"]).tolist()
-  scz_dmn_co = np.ravel(z_values_store["003"]["scz"]["Default-Cingulo-Opercular"]).tolist()
+  con_dmn_co = np.ravel(c_values_store["003"]["con"]["Default-Cingulo-Opercular"]).tolist()
+  scz_dmn_co = np.ravel(c_values_store["003"]["scz"]["Default-Cingulo-Opercular"]).tolist()
 
-  con_fp_cer = np.ravel(z_values_store["003"]["con"]["Fronto-Parietal-Cerebellar"]).tolist()
-  scz_fp_cer = np.ravel(z_values_store["003"]["scz"]["Fronto-Parietal-Cerebellar"]).tolist()
+  con_fp_cer = np.ravel(c_values_store["003"]["con"]["Fronto-Parietal-Cerebellar"]).tolist()
+  scz_fp_cer = np.ravel(c_values_store["003"]["scz"]["Fronto-Parietal-Cerebellar"]).tolist()
 
-  con_dmn_fp = np.ravel(z_values_store["003"]["con"]["Default-Fronto-Parietal"]).tolist()
-  scz_dmn_fp = np.ravel(z_values_store["003"]["scz"]["Default-Fronto-Parietal"]).tolist()
+  con_dmn_fp = np.ravel(c_values_store["003"]["con"]["Default-Fronto-Parietal"]).tolist()
+  scz_dmn_fp = np.ravel(c_values_store["003"]["scz"]["Default-Fronto-Parietal"]).tolist()
 
-  con_fp_co = np.ravel(z_values_store["003"]["con"]["Fronto-Parietal-Cingulo-Opercular"]).tolist()
-  scz_fp_co = np.ravel(z_values_store["003"]["scz"]["Fronto-Parietal-Cingulo-Opercular"]).tolist()
+  con_fp_co = np.ravel(c_values_store["003"]["con"]["Fronto-Parietal-Cingulo-Opercular"]).tolist()
+  scz_fp_co = np.ravel(c_values_store["003"]["scz"]["Fronto-Parietal-Cingulo-Opercular"]).tolist()
 
   # perform permutation test
   dmn_cer_p_value = permute(scz_dmn_cer,con_dmn_cer)  

code/utils/connectivity_utils.py

@@ -6,8 +6,30 @@
 import sys, os
 sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
 import project_config
+import random
 import numpy as np
 
+def permute (r1,r2):
+  """
+  This function performs the permuation test to two lists of r-values (r1:scz; r2:con).
+  Ho: mu_r1 = mu_r2
+  H1: muri < mu_r2
+  input:
+  1.r1 and r2 are two arrays containing the r-values
+  output:
+  1. one sided p-values
+  """
+  n1 = len(r1)
+  n2 = len(r2)
+  t_obs = np.mean(r1)-np.mean(r2)
+  pool = r1+r2
+  diff = []
+  for i in range(0,1000):
+    sample = random.sample(pool,n1)
+    diff.append(np.mean(sample)-(sum(pool)-sum(sample))/n2)
+  p_value = sum(list(i <= t_obs for i in diff))/len(diff)
+  return p_value
+
 def roi_cor (data, roi1,roi2):
   """
   Parameters
@@ -26,7 +48,6 @@ def roi_cor (data, roi1,roi2):
   timecourse2 = [data[roi2[j]] for j in range(0,len(roi2))]
   avg_time2 = np.mean(timecourse2,axis=0)
   cor = np.corrcoef(avg_time1,avg_time2)[1,0]
-
   return cor
 
 def network_cor(data, net1, net2, is_same):

code/utils/tests/test_connectivity_utils.py

@@ -6,7 +6,49 @@
     nosetests connectivity_utils.py
 """
 import numpy as np
+import random
 import os
 from .. import connectivity_utils
 
 from numpy.testing import assert_array_equal, assert_almost_equal
+
+def test_c_within_and_c_between():
+
+  # mocking the correlation values store
+  test_c_values_store = {"test_network_1":{"test_roi_1":((0,0,0),(0,0,1)), "test_roi_2":((0,1,0),(0,1,1))},
+                    "test_network_2":{"test_roi_3":((1,0,0),(1,0,1)), "test_roi_4":((1,1,0),(1,1,1))}}
+
+
+  data = np.zeros((2,2,2,3))
+  data[0,0,0] = [1,2,3]
+  data[0,0,1] = [1,2,3]
+  data[0,1,0] = [-1,-2,-3]
+  data[0,1,1] = [-1,-2,-3]
+  data[1,0,0] = [5,4,37]
+  data[1,0,1] = [5,4,37]
+  data[1,1,0] = [-3,-244,-1]
+  data[1,1,1] = [-3,-244,-1]
+
+
+  actual = connectivity_utils.c_within(data, test_c_values_store)
+
+  # expected values are explicitly calculated according to the rules explained in the paper
+  expected = {'test_network_1':(np.corrcoef([1,2,3],[-1,-2,-3])[1,0],), 'test_network_2': (np.corrcoef([5,4,37],[-3,-244,-1])[1,0],)}
+
+  assert_almost_equal(expected['test_network_1'], expected['test_network_1'])
+  assert_almost_equal(expected['test_network_2'], expected['test_network_2'])
+
+  actual = connectivity_utils.c_between(data, test_c_values_store)
+
+  # expected values are explicitly calculated according to the rules explained in the paper
+  expected = [np.corrcoef([1,2,3],[5,4,37])[1,0], np.corrcoef([1,2,3],[-3,-244,-1])[1,0],np.corrcoef([-1,-2,-3],[5,4,37])[1,0], np.corrcoef([-1,-2,-3],[-3,-244,-1])[1,0]]
+
+  assert_almost_equal(np.sort(expected), np.sort(actual['test_network_1-test_network_2']))
+
+def test_permute():
+  r1 = np.linspace(0,1,num=30).tolist()
+  r2 = np.linspace(0,1,num=20).tolist()
+  random.seed(0)
+  actual = connectivity_utils.permute(r1,r2)
+  expected = 0.496
+  assert_almost_equal(actual,expected)