Structural connectome architecture and biological vulnerability shape cortical atrophy in cocaine use disorder.

1. # %%
2. import os, sys
3. import scipy.io
4. import numpy as np
5. import matplotlib.pyplot as plt
6. from neuromaps.images import load_data
7. from neuromaps.images import dlabel_to_gifti
8. from netneurotools import datasets as nntdata
9. from nilearn.datasets import fetch_atlas_schaefer_2018
10. from statsmodels.stats.multitest import multipletests
11. from scipy.io import loadmat, savemat
12. from neuromaps import images, nulls
15. ref_path = r"D:\OneDrive\References\Farahani_ALS-main\codes" ## https://github.com/netneurolab/Farahani_ALS/tree/main/codes
17. sys.path.append(ref_path)
18. from functions import (parcel2fsLR,
19. save_gifti,
20. vasa_null_Schaefer)
21. from globals import path_results, path_fig, nnodes, path_atlas, path_surface, path_wb_command
23. # %%
24. nspins = 1000 # number of null realizations for the spin test
26. #------------------------------------------------------------------------------
27. # Load group-averegd W-score data of SUDMEX_TMS
28. #------------------------------------------------------------------------------
29. tmp = loadmat(f'D:\\5_TMS_response\\CT_Wscore53avg_S400_7Net_vHCP2.mat')
30. disease_profile = tmp['Wscore_s400_7Net'].squeeze()
32. # %%
33. #------------------------------------------------------------------------------
34. # Load von-economo atlas information
35. #------------------------------------------------------------------------------
36. def load_von_economo_atlas(path_in, nnodes):
37. atlas_data = np.squeeze(scipy.io.loadmat(path_in + 'economo_Schaefer400.mat')['pdata'])
38. return atlas_data - 1, ['primary motor',
39. 'association',
40. 'association',
41. 'primary/secondary sensory',
42. 'primary sensory',
43. 'limbic',
44. 'insular']
45. def load_Yeo7_atlas(path_in, nnodes):
46. atlas_data = np.squeeze(scipy.io.loadmat(path_in + 'Yeo7_Schaefer400.mat')['yeodata'])
47. return atlas_data - 1, ['visual',
48. 'somatomotor',
49. 'dorsal attention',
50. 'salience/ventral attention',
51. 'limbic',
52. 'control',
53. 'default']
57. atlas_7Network_von, label_von_networks = load_von_economo_atlas(path_atlas, nnodes)
58. atlas_7Network_yeo, label_yeo_networks = load_Yeo7_atlas(path_atlas, nnodes)
59. num_labels = len(label_von_networks)
61. # %%
62. import numpy as np
63. from scipy.stats import ttest_1samp
64. from statsmodels.stats.multitest import multipletests
66. def compute_network_ttest_vs_zero(disease_profile, atlas_labels_array, atlas_label_names):
67. """
68. Within each network, a two-tailed one-sample t-test was performed to
69. determine whether its W-score significantly deviated from zero (indicating cortical thinning or thickening)
71. inputs：
72. disease_profile : array, shape (nnodes,)
73. atlas_labels_array : array, shape (nnodes,)
74. atlas_label_names : list of str
77. outputs:
78. result_dict :
79. - "network": list of network names
80. - "mean_value": mean W-score for each network
81. - "p_value": raw p-value from the two-tailed test
82. - "p_fdr": FDR-corrected p-value
83. """
85. num_labels = len(atlas_label_names)
86. mean_values = []
87. p_values = []
89. for i in range(num_labels):
90. # extract the W-score
91. mask = atlas_labels_array == i
92. scores = disease_profile[mask]
93. scores = scores[~np.isnan(scores)] # remove NaN
95. mean_values.append(np.mean(scores))
97. # two-tailed one-sample t-test: H0: μ == 0
98. t_stat, p = ttest_1samp(scores, popmean=0)
99. p_values.append(p)
101. # FDR correction
102. p_fdr = multipletests(p_values, method='fdr_bh')[1]
104. result_dict = {
105. "network": atlas_label_names,
106. "mean_value": mean_values,
107. "p_value": p_values,
108. "p_fdr": p_fdr
109. }
111. return result_dict
113. # %%
114. result = compute_network_ttest_vs_zero(
115. disease_profile=disease_profile,
116. atlas_labels_array=atlas_7Network_yeo,
117. atlas_label_names=label_yeo_networks
118. )
119. for name, mean, p, pfdr in zip(result["network"], result["mean_value"], result["p_value"], result["p_fdr"]):
120. print(f"{name:<20} | Mean = {mean:.4f} | p = {p:.4f} | pFDR = {pfdr:.4f}")
122. # %%
123. result = compute_network_ttest_vs_zero(
124. disease_profile=disease_profile,
125. atlas_labels_array=atlas_7Network_von,
126. atlas_label_names=label_von_networks
127. )
128. for name, mean, p, pfdr in zip(result["network"], result["mean_value"], result["p_value"], result["p_fdr"]):
129. print(f"{name:<20} | Mean = {mean:.4f} | p = {p:.4f} | pFDR = {pfdr:.4f}")
131. # %%
132. ##visualization the results
133. def plot_boxplot(nulls,path_fig=None, name_to_save=None):
134. """
135. Create a box plot for the null distribution
136. and visualize the actual value on top
137. """
138. fig, axes = plt.subplots(figsize = (7, 4))
139. axes.boxplot(nulls, vert=True,flierprops=dict(markersize=6))
141. axes.spines['top'].set_visible(False)
142. axes.spines['right'].set_visible(False)
143. axes.spines['bottom'].set_visible(False)
145. plt.tight_layout()
146. if name_to_save:
147. plt.savefig(os.path.join(path_fig, name_to_save + '.pdf'),dpi=300)
148. plt.show()
151. # Get actual/real values for each network YEO
152. network_specific_disease_measure_yeo = []
153. for label_ind in range(num_labels):
154. temp = []
155. for roi_ind in range(nnodes):
156. if atlas_7Network_yeo[roi_ind] == label_ind:
157. temp.append(disease_profile[roi_ind])
158. print(np.array(temp).shape)
159. network_specific_disease_measure_yeo.append(np.array(temp))
161. # %%
162. ## sort the data YEO
163. row_medians = [np.nanmedian(arr) for arr in network_specific_disease_measure_yeo]
164. sorted_indices_Yeo = np.argsort(row_medians)[::-1]
165. sorted_Yeo_data_TMS = [network_specific_disease_measure_yeo[i] for i in sorted_indices_Yeo]
166. plot_boxplot(sorted_Yeo_data_TMS,path_fig='D:\\OneDrive\\5_TMS_response\\',
167. name_to_save='CT_Wscore53avg_S400_7Net_vHCP2_Yeo_hor')
169. # %%
170. # Get actual/real values for each network VON
171. network_specific_disease_measure_von = []
172. for label_ind in range(num_labels):
173. temp = []
174. for roi_ind in range(nnodes):
175. if atlas_7Network_von[roi_ind] == label_ind:
176. temp.append(disease_profile[roi_ind])
177. print(np.array(temp).shape)
178. network_specific_disease_measure_von.append(np.array(temp))
180. # %%
181. ## sort the data VON
182. row_medians = [np.nanmedian(arr) for arr in network_specific_disease_measure_von]
183. sorted_indices_von = np.argsort(row_medians)[::-1]
184. sorted_von_data_TMS = [network_specific_disease_measure_von[i] for i in sorted_indices_von]
185. plot_boxplot(sorted_von_data_TMS,path_fig='D:\\OneDrive\\5_TMS_response\\',
186. name_to_save='CT_Wscore53avg_S400_7Net_vHCP2_von_hor')

S1_Network_atrophy_stats.ipynb at commit 9e32f25, no license · at the source