OSCR

Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning.

Code ↔ Paper

7 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 7 matches
  1. [1] § Results › GMV ROI ~ Behaviour Correlation Analysis ↔ ROI_Bx_T1_correlation.ipynb, lines 744–783 · score 0.66 · GMV change scores, Block12 Block3, d5 d2, LOSO, CI, Pearsons
  2. [2] § Statistical Analyses › ROI Analyses ↔ ROI_Bx_T1_correlation.ipynb, lines 744–783 · score 0.65 · Block12 Block3, change scores, GMV change, Pearson, correlation, behavioural
  3. [3] § Statistical Analyses › ROI Analyses ↔ ROI_Bx_T1_correlation.ipynb, lines 457–567 · score 0.60 · GM mask, GM ROI, WM, TFCE, correlation, clusters
  4. [4] § Methods › Voxel‐Based Morphometry ↔ longitudinal_segmentation_all40subjects.m, lines 14–49 · score 0.58 · shooting, template, CAT12, modulated, SPM12, segmentation
  5. [5] § Statistical Analyses › GM Volume Change During MSL Across Learning Stages ↔ ROI_Bx_T1_correlation.ipynb, lines 457–567 · score 0.57 · GM voxels, modulated, threshold, mask, tissue, TFCE
  6. [6] § Results › T1 ROI ~ Behaviour Relationship ↔ ROI_Bx_T1_correlation.ipynb, lines 605–699 · score 0.55 · SYN changes, left SPC, modalities, CI, VBM, correlation
  7. [7] § Statistical Analyses › GM Volume Change During MSL Across Learning Stages ↔ flexible_factorial_design_withbinmask.m, lines 421–507 · score 0.53 · flexible factorial design, threshold, mask, voxels, d1, d2

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 875 lines · 37 KB · no license · 5 matches

  1. # %%
  2. import os
  3. from os.path import join
  4. from glob import glob
  5. import csv
  7. import nibabel as nb
  8. import numpy as np
  9. import nilearn
  10. import pandas as pd
  12. from nilearn.image import threshold_img
  13. from nilearn import image, plotting
  14. from nilearn.maskers import NiftiMasker
  16. from scipy.ndimage import binary_dilation
  17. from scipy import stats
  18. from scipy.stats import pearsonr,chi2
  19. import scipy.stats.mstats as mstats
  21. import statsmodels.formula.api as smf
  22. from statsmodels.stats.multitest import multipletests
  23. import statsmodels.api as sm
  24. from statsmodels.stats.anova import AnovaRM
  25. from statsmodels.stats.multitest import multipletests
  27. from matplotlib import pyplot as plt
  28. import seaborn as sns
  30. import warnings
  31. warnings.filterwarnings('ignore')
  32. from sklearn.model_selection import LeaveOneOut
  34. # %%
  35. ##### Fisher's z test for T1 and Bx between LRN adn SMP correlations.#####
  36. ##### LRN n=19 and SMP n= 20
  37. z1 = 0.5 * np.log((1 + -0.65) / (1 - -0.65))
  38. z2 = 0.5 * np.log((1 + 0.18) / (1 - 0.18))
  39. SE = np.sqrt(1/(19-3) + 1/(20-3))
  40. Z = (z1 - z2) / SE
  41. p = 2 * (1 - stats.norm.cdf(abs(Z)))
  43. print(f"Z = {Z:.3f}, p = {p:.4f}")
  45. # %%
  46. #### Effect sizes for Table 2 of ms ######
  47. data1 = {
  48. 'LRN_mean': [0.0622, 0.0470, 0.0639, 0.1181, 0.0690],
  49. 'LRN_sem': [0.0210, 0.0375, 0.0244, 0.0631, 0.0455],
  50. 'SMP_mean': [-0.0558, -0.1071, -0.0819, -0.0890, -0.1097],
  51. 'SMP_sem': [0.0233, 0.0348, 0.0195, 0.0312, 0.0269]
  52. }
  54. df1 = pd.DataFrame(data1)
  55. n_lrn = 19
  56. n_smp = 20
  58. # Convert SEM to SD
  59. df1['LRN_sd'] = df1['LRN_sem'] * np.sqrt(n_lrn)
  60. df1['SMP_sd'] = df1['SMP_sem'] * np.sqrt(n_smp)
  62. # Weighted pooled SD (correct for unequal n)
  63. df1['SD_pooled'] = np.sqrt(
  64. (((n_lrn - 1) * df1['LRN_sd']**2) + ((n_smp - 1) * df1['SMP_sd']**2)) /
  65. (n_lrn + n_smp - 2)
  66. )
  68. # Mean difference
  69. df1['Mean_diff'] = df1['LRN_mean'] - df1['SMP_mean']
  71. # Cohen's d
  72. df1['Cohens_d'] = df1['Mean_diff'] / df1['SD_pooled']
  74. # Display
  75. print(df1[['LRN_mean', 'SMP_mean', 'Mean_diff', 'Cohens_d']].round(2))
  77. # %%
  78. df = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/SYN_Scores_by_block_Jhelum.csv')
  79. df = df[df['PID'] != 'P11']
  81. n_lrn = df[df['Group'] == 'LRN']['PID'].nunique()
  82. n_smp = df[df['Group'] == 'SMP']['PID'].nunique()
  84. day_columns = ['Day1', 'Day2', 'Day3', 'Day4', 'Day5', 'Day17']
  85. lrn_averages = []
  86. smp_averages = []
  87. sem_LRN =[]
  88. sem_SMP =[]
  90. for day in day_columns:
  91. lrn_avg = df[df['Group'] == 'LRN'][day].mean()
  92. lrn_checkstd = df[df['Group'] == 'LRN'][day].std()
  93. lrn_checksem = lrn_checkstd/np.sqrt(n_lrn-1)
  94. smp_avg = df[df['Group'] == 'SMP'][day].mean()
  95. smp_checkstd = df[df['Group'] == 'SMP'][day].std()
  96. smp_checksem = smp_checkstd/np.sqrt(n_smp-1)
  97. lrn_averages.append(lrn_avg)
  98. smp_averages.append(smp_avg)
  99. sem_LRN.append(lrn_checksem)
  100. sem_SMP.append(smp_checksem)
  102. plt.figure(figsize=(10,7))
  103. day_split_index = day_columns.index('Day5')
  104. plt.errorbar(day_columns[:day_split_index + 1],
  105. lrn_averages[:day_split_index + 1],
  106. yerr=sem_LRN[:day_split_index + 1],
  107. capsize=3, marker='o', label='LRN', linewidth=2, color='#004586')
  108. plt.errorbar(day_columns[:day_split_index + 1],
  109. smp_averages[:day_split_index + 1],
  110. yerr=sem_SMP[:day_split_index + 1],
  111. capsize=3, marker='s', label='SMP', linewidth=2, color='#FF420E')
  112. plt.errorbar(day_columns[day_split_index:],
  113. lrn_averages[day_split_index:],
  114. yerr=sem_LRN[day_split_index:],
  115. capsize=3, marker='o', linestyle='dotted', linewidth=2, color='#004586')
  116. plt.errorbar(day_columns[day_split_index:],
  117. smp_averages[day_split_index:],
  118. yerr=sem_SMP[day_split_index:],
  119. capsize=3, marker='s', linestyle='dotted', linewidth=2, color='#FF420E')
  120. plt.legend(fontsize=14)
  121. plt.xlabel('Day', fontsize=18)
  122. plt.ylabel('SYN Score (in ms)', fontsize=18)
  123. plt.ylim(0,230)
  124. plt.yticks(np.arange(0, 201, 50), fontsize=16)
  125. plt.xticks(fontsize=16)
  126. plt.tight_layout()
  127. #plt.savefig('bx_adjusted_current.png', dpi=600)
  128. plt.show()
  130. # %%
  132. # df is already loaded with dataset for SYN block by block, excluding the outlier subject of LRN group P11
  134. data_long = df.melt(
  135. id_vars=['PID', 'Group'],
  136. value_vars=day_columns,
  137. var_name='Day',
  138. value_name='SYN')
  139. data_long['SYN'] = pd.to_numeric(data_long['SYN'], errors='coerce')
  140. data_long = data_long.dropna(subset=['SYN'])
  141. data_long['Group'] = data_long['Group'].astype('category')
  142. data_long['Day'] = pd.Categorical(data_long['Day'],
  143. categories=['Day1','Day2','Day3','Day4','Day5','Day17'],
  144. ordered=True)
  145. data_long['PID'] = data_long['PID'].astype('category')
  146. print(f"\nSample sizes:")
  147. print(f" LRN: {data_long[data_long['Group']=='LRN']['PID'].nunique()} participants")
  148. print(f" SMP: {data_long[data_long['Group']=='SMP']['PID'].nunique()} participants")
  149. print(f" Total: {data_long['PID'].nunique()} participants")
  150. print(f" Total rows: {len(data_long)}\n")
  152. lrn = data_long[data_long['Group'] == 'LRN']
  153. pairs = [('Day1','Day2'), ('Day2','Day3'), ('Day3','Day4'), ('Day4','Day5'), ('Day5','Day17')]
  155. pvals, tvals, gvals = [], [], []
  157. def hedges_g(x, y):
  158. n = len(x)
  159. d = (np.mean(x) - np.mean(y)) / np.std(np.concatenate([x, y]), ddof=1)
  160. return d * (1 - (3 / (4 * n - 9)))
  162. for a, b in pairs:
  163. vals_a = lrn[lrn['Day'] == a]['SYN'].values
  164. vals_b = lrn[lrn['Day'] == b]['SYN'].values
  165. t, p = stats.ttest_rel(vals_a, vals_b)
  166. g = hedges_g(vals_a, vals_b)
  167. tvals.append(t)
  168. pvals.append(p)
  169. gvals.append(g)
  171. rej, p_corr, _, _ = multipletests(pvals, alpha=0.05, method='bonferroni')
  173. for i, (a, b) in enumerate(pairs):
  174. sig = '***' if p_corr[i] < 0.001 else '**' if p_corr[i] < 0.01 else '*' if p_corr[i] < 0.05 else 'n.s.'
  175. print(f"{a} vs {b}: t={tvals[i]:6.2f}, p_raw={pvals[i]:.4f}, p_corr={p_corr[i]:.4f} ({sig}), Hedges' g={gvals[i]:5.2f}")
  178. # %%
  179. # COMPREHENSIVE DESCRIPTIVE STATISTICS
  180. # Calculate statistics for each Group × Day combination
  181. # Filter the data_long to exclude the 'LRN SMP task' group
  182. data_long_filtered = data_long[data_long['Group'].isin(['LRN', 'SMP'])].copy()
  183. data_long_filtered['Group'] = data_long_filtered['Group'].cat.remove_unused_categories()
  184. stats_summary = []
  185. for group in ['LRN', 'SMP']:
  186. for day in ['Day1', 'Day2', 'Day3', 'Day4', 'Day5', 'Day17']:
  187. subset = data_long_filtered[(data_long_filtered['Group'] == group) &
  188. (data_long_filtered['Day'] == day)]['SYN']
  190. if len(subset) > 0:
  191. stats_dict = {
  192. 'Group': group,
  193. 'Day': day,
  194. 'N': len(subset),
  195. 'Mean': subset.mean(),
  196. 'SD': subset.std(),
  197. 'SEM': subset.sem(),
  198. 'Median': subset.median(),
  199. 'Q1': subset.quantile(0.25),
  200. 'Q3': subset.quantile(0.75),
  201. 'IQR': subset.quantile(0.75) - subset.quantile(0.25),
  202. 'Min': subset.min(),
  203. 'Max': subset.max(),
  204. 'Range': subset.max() - subset.min(),
  205. 'Skewness': stats.skew(subset),
  206. 'Kurtosis': stats.kurtosis(subset)
  207. }
  208. stats_summary.append(stats_dict)
  210. stats_df = pd.DataFrame(stats_summary)
  211. print("\n CENTRAL TENDENCY AND VARIABILITY")
  212. for group in ['LRN', 'SMP']:
  213. print(f"\n{group} GROUP:")
  214. group_data = stats_df[stats_df['Group'] == group]
  215. print(f"{'Day':<8} {'N':<4} {'Mean±SD':<20} {'Median (IQR)':<25} {'Range':<15}")
  216. print("-" * 72)
  217. for _, row in group_data.iterrows():
  218. mean_sd = f"{row['Mean']:.2f} ± {row['SD']:.2f}"
  219. median_iqr = f"{row['Median']:.2f} ({row['Q1']:.2f}-{row['Q3']:.2f})"
  220. range_str = f"{row['Min']:.2f}-{row['Max']:.2f}"
  221. print(f"{row['Day']:<8} {int(row['N']):<4} {mean_sd:<20} {median_iqr:<25} {range_str:<15}")
  222. print("\n DISTRIBUTION CHARACTERISTICS")
  223. for group in ['LRN', 'SMP']:
  224. print(f"\n{group} GROUP:")
  225. group_data = stats_df[stats_df['Group'] == group]
  226. print(f"{'Day':<8} {'Skewness':<12} {'Kurtosis':<12} {'Distribution Shape':<30}")
  227. print("-" * 62)
  228. for _, row in group_data.iterrows():
  229. if abs(row['Skewness']) < 0.5:
  230. skew_interp = "Approximately symmetric"
  231. elif row['Skewness'] > 0:
  232. skew_interp = "Right-skewed (positive)"
  233. else:
  234. skew_interp = "Left-skewed (negative)"
  235. if abs(row['Kurtosis']) < 0.5:
  236. kurt_interp = "Mesokurtic (normal)"
  237. elif row['Kurtosis'] > 0:
  238. kurt_interp = "Leptokurtic (heavy-tailed)"
  239. else:
  240. kurt_interp = "Platykurtic (light-tailed)"
  242. distribution = f"{skew_interp}, {kurt_interp}"
  243. print(f"{row['Day']:<8} {row['Skewness']:>6.2f} {row['Kurtosis']:>6.2f} {distribution}")
  246. print("\n OVERALL SUMMARY ACROSS ALL DAYS")
  247. for group in ['LRN', 'SMP']:
  248. group_all = data_long_filtered[data_long_filtered['Group'] == group]['SYN']
  249. print(f"\n{group} GROUP (N={len(group_all)} observations):")
  250. print(f" Overall Mean ± SD: {group_all.mean():.2f} ± {group_all.std():.2f} ms")
  251. print(f" Overall Median: {group_all.median():.2f} ms")
  252. print(f" Overall Range: {group_all.min():.2f} - {group_all.max():.2f} ms")
  253. print(f" Interquartile Range: {group_all.quantile(0.25):.2f} - {group_all.quantile(0.75):.2f} ms")
  255. palette = {
  256. 'LRN': '#004586',
  257. 'SMP': '#FF420E'
  258. }
  260. sns.set(style='whitegrid', font_scale=1.3)
  261. plt.figure(figsize=(14, 7))
  264. sns.violinplot(data=data_long_filtered, x='Day', y='SYN', hue='Group', palette=palette,
  265. inner=None, cut=0, alpha=0.5)
  266. sns.boxplot(data=data_long_filtered, x='Day', y='SYN', hue='Group', palette=palette,
  267. showcaps=False, boxprops={'facecolor':'none', 'zorder':2},
  268. showfliers=False, whiskerprops={'linewidth':0})
  269. sns.stripplot(data=data_long_filtered, x='Day', y='SYN', hue='Group',
  270. dodge=True, jitter=True, alpha=0.5, size=5, color='k')
  271. plt.title('Raincloud Plot of SYN Scores - LRN vs SMP', fontsize=16)
  272. plt.xlabel('Training Day', fontsize=16)
  273. plt.ylabel('SYN (ms)', fontsize=16)
  274. handles, labels = plt.gca().get_legend_handles_labels()
  275. plt.legend(handles[:2], labels[:2], title='Group', bbox_to_anchor=(1.05, 1), loc='upper left')
  276. plt.tight_layout()
  277. #plt.savefig('raincloud_graph_allsubs_current.png', dpi=600, bbox_inches='tight')
  278. plt.show()
  281. # %% [markdown]
  282. # # DF FOR VBM CLUSTERS
  284. # %%
  285. #cluster D- slow learning L SPC
  286. dfvbm_LRN_D = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/clusterD_VBM.csv',usecols=['LRN_d1','LRN_d2','LRN_d5', 'LRN_d6'])
  287. dfvbm_SMP_D = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/clusterD_VBM.csv',usecols=['SMP_d1','SMP_d2','SMP_d5','SMP_d6'])
  289. synscores_LRN = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/syn_score_25feb2023.csv',usecols=['Block1','Block3','Block12'])
  291. synB12minusB3 = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/syn_score_25feb2023.csv',usecols=['B12minusB3']) #d2 vs d5 slow learning -> Cluster D LEFT SPC
  293. synB3minusB1 = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/syn_score_25feb2023.csv',usecols=['B3minusB1']) #d2 vs d5 slow learning -> Cluster D LEFT SPC
  295. synscores_SMP = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/SYN_BLOCKS_SMP.csv',usecols=['Block 1','Block 3','Block 12','B12minusB3'])
  297. # %%
  298. # Change in vbm for learning stages
  299. dfvbm_LRN_D['d5minusd2'] = dfvbm_LRN_D['LRN_d5'] - dfvbm_LRN_D['LRN_d2'] #cluster D- slow learning
  300. new_dfvbm_LRN_D = dfvbm_LRN_D['d5minusd2'] #cluster D- slow learning
  301. var_vec_d = new_dfvbm_LRN_D.copy()
  302. var_vec_d = var_vec_d[~np.isnan(var_vec_d)]
  303. pearsonr( var_vec_d, synB12minusB3['B12minusB3'])
  305. # %%
  306. # SMP group GMV vs SYN
  307. dfvbm_SMP_D['d5minusd2'] = dfvbm_SMP_D['SMP_d5'] - dfvbm_SMP_D['SMP_d2']
  308. new_dfvbm_SMP_D = dfvbm_SMP_D['d5minusd2'] #cluster D- slow learning
  309. var_vec_d_SMP = new_dfvbm_SMP_D.copy()
  310. var_vec_d_SMP = var_vec_d_SMP[~np.isnan(var_vec_d_SMP)]
  311. pearsonr( var_vec_d_SMP, synscores_SMP['B12minusB3'])
  313. # %%
  314. slope2, intercept2, r2, p2, stderr2 = mstats.linregress(var_vec_d, synB12minusB3['B12minusB3'])
  315. line2 = f'Regression line: y={intercept2:.2f}+{slope2:.2f}x'
  317. r2_value = r2
  318. p2_value = p2
  320. # Calculate CI
  321. def correlation_ci(r, n, confidence=0.95):
  322. """Calculate confidence interval for correlation coefficient"""
  323. from scipy.stats import norm
  324. z = np.arctanh(r)
  325. se = 1 / np.sqrt(n - 3)
  326. z_crit = norm.ppf((1 + confidence) / 2)
  327. ci_lower_z = z - z_crit * se
  328. ci_upper_z = z + z_crit * se
  329. ci_lower = np.tanh(ci_lower_z)
  330. ci_upper = np.tanh(ci_upper_z)
  331. return ci_lower, ci_upper
  333. n_samples = len(var_vec_d)
  334. ci_lower, ci_upper = correlation_ci(r2_value, n_samples)
  336. print(f"r = {r2_value:.2f}, 95% CI [{ci_lower:.2f}, {ci_upper:.2f}], p = {p2_value:.3f}, n = {n_samples}")
  338. line2 = f'Regression line: y={intercept2:.2f}+{slope2:.2f}x'
  339. r2 = f'r = {r2_value:.2f}'
  340. p2 = f'p = {p2_value:.2f}'
  342. # %%
  343. x_vals2 = np.array([var_vec_d.min(), var_vec_d.max()])
  344. y_vals2 = intercept2 + slope2 * x_vals2
  345. angle2 = np.degrees(np.arctan2(y_vals2[1] - y_vals2[0], x_vals2[1] - x_vals2[0]))
  346. x_text2 = x_vals2[0] + 0.054 # adjust as needed
  347. y_text2 = (intercept2+ 8) + slope2 * x_text2
  348. x_text3 = x_vals2[0] + 0.055 # adjust as needed
  349. y_text3 = (intercept2+ -15) + slope2 * x_text3
  351. fig, ax = plt.subplots()
  352. ax.plot(var_vec_d, synB12minusB3['B12minusB3'], linewidth=0, marker='o', label='Participants')
  353. ax.plot(var_vec_d, intercept2 + slope2 * var_vec_d, label=line2)
  354. ax.set_xlabel('Change in mean GM (d5-d2)',fontsize=18)
  355. ax.set_ylabel('SYN score (B12-B3)',fontsize=18)
  356. ax.legend(loc= 2, facecolor='white')
  357. plt.title('Left Superior Parietal Cortex', fontsize=18)
  358. plt.text(x_text2, y_text2, r2, fontsize=12, rotation=angle2, rotation_mode='anchor', transform_rotates_text=True)
  359. plt.text(x_text3, y_text3, p2, fontsize=12, rotation=angle2, rotation_mode='anchor', transform_rotates_text=True)
  360. ax.yaxis.grid(True) # Enables horizontal grid lines
  361. ax.xaxis.grid(False)
  362. plt.ylim(-150,100)
  363. plt.xticks(fontsize=14)
  364. plt.yticks(fontsize=16)
  365. plt.tight_layout()
  366. #plt.savefig('L_SPC_slow_cluster_D.png', dpi=600)
  367. plt.show()
  369. # %% [markdown]
  370. # # T1 mean value extraction from all subjects and each day
  372. # %%
  373. ## load the subject ids here ##
  374. LRN = ['P05','P06','P07','P08','P09', 'P10', 'P12', 'P14', 'P16', 'P17','P18', 'P20', 'P21', 'P26', 'P28','P31', 'P36', 'P37', 'P38']
  375. SMP = ['P13','P15','P22','P23','P24', 'P25', 'P27', 'P30', 'P32', 'P33', 'P34','P35', 'P39', 'P40', 'P41', 'P42','P43', 'P44', 'P45', 'P46']
  377. # %%
  378. #Gather the participants file names for GM and WM d1/d2/d5 in LRN group
  380. Pname_LRN_GM_d1 = sorted(list(set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P*/d1/mri/mwp1r*.nii"))-set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P11/d1/mri/mwp1r*.nii"))))
  381. Pname_LRN_GM_d2 = sorted(list(set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P*/d2/mri/mwp1r*.nii"))-set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P11/d2/mri/mwp1r*.nii"))))
  382. Pname_LRN_GM_d5 = sorted(list(set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P*/d5/mri/mwp1r*.nii"))-set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P11/d5/mri/mwp1r*.nii"))))
  383. Pname_LRN_WM_d1 = sorted(list(set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P*/d1/mri/mwp2r*.nii"))-set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P11/d1/mri/mwp2r*.nii"))))
  384. Pname_LRN_WM_d2 = sorted(list(set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P*/d2/mri/mwp2r*.nii"))-set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P11/d2/mri/mwp2r*.nii"))))
  385. Pname_LRN_WM_d5 = sorted(list(set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P*/d5/mri/mwp2r*.nii"))-set(glob("/data/neuralabc/paujhe/mMPI/data/Learning/P11/d5/mri/mwp2r*.nii"))))
  387. #Gather the participants file names for GM and WM d1/d2/d5 in SMP group
  389. Pname_SMP_GM_d1 = sorted(glob("/data/neuralabc/paujhe/mMPI/data/Control/P*/d1/mri/mwp1r*.nii"))
  390. Pname_SMP_GM_d2 = sorted(glob("/data/neuralabc/paujhe/mMPI/data/Control/P*/d2/mri/mwp1r*.nii"))
  391. Pname_SMP_GM_d5 = sorted(glob("/data/neuralabc/paujhe/mMPI/data/Control/P*/d5/mri/mwp1r*.nii"))
  392. Pname_SMP_WM_d1 = sorted(glob("/data/neuralabc/paujhe/mMPI/data/Control/P*/d1/mri/mwp2r*.nii"))
  393. Pname_SMP_WM_d2 = sorted(glob("/data/neuralabc/paujhe/mMPI/data/Control/P*/d2/mri/mwp2r*.nii"))
  394. Pname_SMP_WM_d5 = sorted(glob("/data/neuralabc/paujhe/mMPI/data/Control/P*/d5/mri/mwp2r*.nii"))
  397. # %%
  398. #Gather the participants file names for FA/MD d1/d2/d5 in LRN group
  399. files1 = []
  400. files2 = []
  401. files3 = []
  402. files4 = []
  403. files5 = []
  404. files6 = []
  405. for participant in LRN:
  406. files1.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/FA/{participant}_d1_*_FA.nii"))
  407. files2.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/FA/{participant}_d2_*_FA.nii"))
  408. files3.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/FA/{participant}_d5_*_FA.nii"))
  409. files4.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/MD/{participant}_d1_*_MD.nii"))
  410. files5.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/MD/{participant}_d2_*_MD.nii"))
  411. files6.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/MD/{participant}_d5_*_MD.nii"))
  413. Pname_LRN_FA_d1 = sorted(files1)
  414. Pname_LRN_FA_d2 = sorted(files2)
  415. Pname_LRN_FA_d5 = sorted(files3)
  416. Pname_LRN_MD_d1 = sorted(files4)
  417. Pname_LRN_MD_d2 = sorted(files5)
  418. Pname_LRN_MD_d5 = sorted(files6)
  420. # %%
  421. #Gather the participants file names for FA/MD d1/d2/d5 in SMP group
  422. files1 = []
  423. files2 = []
  424. files3 = []
  425. files4 = []
  426. files5 = []
  427. files6 = []
  428. for participant in SMP:
  429. files1.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/FA/{participant}_d1_*_FA.nii"))
  430. files2.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/FA/{participant}_d2_*_FA.nii"))
  431. files3.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/FA/{participant}_d5_*_FA.nii"))
  432. files4.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/MD/{participant}_d1_*_MD.nii"))
  433. files5.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/MD/{participant}_d2_*_MD.nii"))
  434. files6.extend(glob(f"/data/neuralabc/treste/mMPI_DWI/processing/SPM_Analyses/All_maps_MNI/MD/{participant}_d5_*_MD.nii"))
  436. Pname_SMP_FA_d1 = sorted(files1)
  437. Pname_SMP_FA_d2 = sorted(files2)
  438. Pname_SMP_FA_d5 = sorted(files3)
  439. Pname_SMP_MD_d1 = sorted(files4)
  440. Pname_SMP_MD_d2 = sorted(files5)
  441. Pname_SMP_MD_d5 = sorted(files6)
  443. # %%
  444. #Gather the participants file names for qT1map d1/d2/d5 in LRN group
  446. Pname_LRN_T1_d1 = sorted(list(set(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR1_*_d1_*.nii"))-set(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR1_P11_d1_*.nii"))))
  447. Pname_LRN_T1_d2 = sorted(list(set(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR1_*_d2_*.nii"))-set(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR1_P11_d2_*.nii"))))
  448. Pname_LRN_T1_d5 = sorted(list(set(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR1_*_d5_*.nii"))-set(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR1_P11_d5_*.nii"))))
  450. #Gather the participants file names for qT1map d1/d2/d5 in SMP group
  452. #Gather the participants file names for qT1map d1/d2/d5 in LRN group
  453. Pname_SMP_T1_d1 = sorted(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR2_*_d1_*.nii"))
  454. Pname_SMP_T1_d2 = sorted(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR2_*_d2_*.nii"))
  455. Pname_SMP_T1_d5 = sorted(glob("/data/neuralabc/giachi/mMPI/processing/T1m/brain/GR2_*_d5_*.nii"))
  457. # %%
  458. def metric_extract_ROI_gm_wm_gmwm_cut(metric_f, original_ROI_f, gm_seg_f, wm_seg_f, tfce_sig_results,tfce_cut=1.3,
  459. orig_ROI_dil_iterations=3,ext_WM_dil_iterations=1, cut=0,verbosity=0,working_dir='./',
  460. return_testing_images=False):
  461. '''
  462. ext_WM_dil_iterations should likely stay between 1 and 2, CHECK YOUR OUTPUT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! <---------------------------!!!!
  463. Both dil_iterations must be >0
  464. ext_WM_dil_iterations will remove voxels from gm_wm_interface if grows too large (as we are not using a distance based approach)
  466. '''
  468. # the input GM mask is actually a modulated VBM, so we take it at >0 to form a binary GM mask, and use this to then define WM
  469. import subprocess
  470. from os.path import join, basename
  472. #convert all masks to the same space as the metric, store in working_dir
  473. tmp_ROI_fname = join(working_dir,'XXX_tmp_ROI.nii.gz')
  474. cmd = ['mrgrid',original_ROI_f, 'regrid','-template',metric_f,'-interp','nearest',tmp_ROI_fname,'-force']
  475. subprocess.run(cmd)
  477. tmp_gm_fname = join(working_dir,'XXX_tmp_gm.nii.gz')
  478. cmd = ['mrgrid',gm_seg_f, 'regrid','-template',metric_f,'-interp','nearest',tmp_gm_fname,'-force']
  479. subprocess.run(cmd)
  481. tmp_wm_fname = join(working_dir,'XXX_tmp_wm.nii.gz')
  482. cmd = ['mrgrid',wm_seg_f, 'regrid','-template',metric_f,'-interp','nearest',tmp_wm_fname,'-force']
  483. subprocess.run(cmd)
  485. tmp_tfce_fname = join(working_dir,'XXX_tmp_tfce.nii.gz')
  486. cmd = ['mrgrid',tfce_sig_results, 'regrid','-template',metric_f,'-interp','nearest',tmp_tfce_fname,'-force']
  487. subprocess.run(cmd)
  489. # load masks and ROI
  490. ROI_img = nb.load(tmp_ROI_fname)
  491. ROI = ROI_img.get_fdata().astype(bool)
  492. ROI_dil = binary_dilation(ROI,iterations=orig_ROI_dil_iterations) #iteration is what we include when calling the function
  495. #tfce sig results, convert to mask
  496. tfce_mask = nb.load(tmp_tfce_fname).get_fdata()>tfce_cut #tfce_cut is fixed at 1.3 to get p<0.05
  498. # load metric data
  499. metric_d = nb.load(metric_f).get_fdata()
  501. gm = nb.load(tmp_gm_fname).get_fdata()
  502. gm = gm>cut # taking all the values in the gm which are above cut, eg.cut here is 0.
  504. gm_ROI = (gm)*ROI_dil*tfce_mask #these are the voxels in sig GM
  505. gm_ROI_dil = binary_dilation((gm_ROI),iterations=1)
  507. # gm_dil = binary_dilation((gm*ROI_dil),iterations=1) #this dilation count is different to the ROI dilation, only used for gm_wm_int and wm
  509. wm = nb.load(tmp_wm_fname).get_fdata()
  510. wm = wm>0 #this is very liberal if using vbm outputs to be sure everything is captured.
  512. # wm = binary_dilation((wm*ROI_dil),iterations=1) #JP
  513. # gm_wm_int = gm_dil*np.logical_not(gm) * np.logical_not(wm) #JP
  514. gm_wm_int = gm_ROI_dil * np.logical_not(gm_ROI) * wm
  515. wm = wm * np.logical_not(gm_wm_int)*np.logical_not(gm) #refine WM to remove interface
  516. gm_wm_ROI = gm_wm_int # this is the interface, NOT masked by tfce region
  517. # gm_wm_ROI = gm_wm_int*ROI_dil*tfce_mask #interface mask within tfce region
  520. # wm_ROI_tfce = (wm)*tfce_mask*np.logical_not(gm_wm_ROI) #mask by tfce, likely to have very few voxels (maybe none)
  521. wm_ROI_ext = (wm)*binary_dilation((gm_wm_ROI),iterations=ext_WM_dil_iterations)*np.logical_not(gm_wm_ROI) #wm external to the GM, no matter if it is in the tfce or not
  522. gm_wm_ROI = gm_wm_ROI * (binary_dilation(wm_ROI_ext,iterations=1)*np.logical_not(wm_ROI_ext)) #refine based on external WM to remove any potential GM voxels b/c of tissue thresholding issues
  524. ## WM_ROI_ext extraction is based on double dilated gm_ROI
  525. ROI_dil_subset = ROI_dil*tfce_mask #this is just extracting ROI_dil within the tfce significant cluster
  527. # gwm_ROI = (((gm>0)*(~(gm_ROI)))*((wm>0)*(~(wm_ROI))))*ROI_dil #THIS INTERFACE IS NOT CORRECT, can include other regions as well
  529. if verbosity>0:
  530. print(ROI_dil.sum())
  531. print('number vox above cutoff in ROI for gm and wm')
  532. print(np.sum(gm[ROI_dil]))
  533. print(np.sum(wm[ROI_dil]))
  534. print(np.sum(gm_wm_int[ROI_dil]))
  536. print('after masking by sig result')
  537. print(f'gm:\t{gm_ROI.sum()}')
  538. print(f'wm_ext:\t{wm_ROI_ext.sum()}')
  539. # print(f'wm_tfce:\t{wm_ROI_tfce.sum()}')
  540. print(f'gm_wm:\t{gm_wm_ROI.sum()}')
  541. # print(gwm_ROI.sum())
  543. gm_cnt = (gm_ROI.sum())
  544. # wm_cnt = (wm_ROI_tfce.sum())
  545. wm_ext_cnt = (wm_ROI_ext.sum())
  546. gm_wm_cnt = (gm_wm_ROI.sum())
  548. #do they share the same vox?
  549. #this can never happen in this implementation, as we explicitly remove GM from WM
  550. if ((gm_ROI*wm_ROI_ext).sum() > 0):
  551. print("There is overlap between WM and GM, this should not happen!!!!!")
  552. return None
  553. elif ((wm_ROI_ext*gm_wm_ROI).sum() > 0) or ((gm_ROI*gm_wm_ROI).sum() > 0):
  554. print("There is overlap between WM/GM interface and WM and/or GM, this should not happen!!!!!")
  555. return None
  556. print(tmp_ROI_fname)
  557. img_gm = nb.Nifti1Image(gm_ROI,affine=ROI_img.affine,header=ROI_img.header)
  558. # img_wm = nb.Nifti1Image(wm_ROI_tfce,affine=ROI_img.affine,header=ROI_img.header) #thresholded by tfce sig result
  559. img_wm_ext = nb.Nifti1Image(wm_ROI_ext,affine=ROI_img.affine,header=ROI_img.header) #defined only in WM, not thresholded by tfce sig result
  560. img_gmwm = nb.Nifti1Image(gm_wm_ROI,affine=ROI_img.affine,header=ROI_img.header)
  561. if return_testing_images:
  562. return img_gm,img_wm_ext, img_gmwm
  563. else:
  564. return {'ROI_dil_cnt':ROI_dil_subset.sum(), 'ROI_dil_mean':np.mean(metric_d[ROI_dil_subset]),
  565. 'gm_cnt':gm_cnt, 'gm_mean':np.mean(metric_d[gm_ROI]),
  566. 'wm_ext_cnt':wm_ext_cnt, 'wm_ext_mean':np.mean(metric_d[wm_ROI_ext]),
  567. 'gm_wm_interface_cnt':gm_wm_cnt, 'gm_wm_interface_mean':np.mean(metric_d[gm_wm_ROI])}
  569. # %%
  570. original_ROI_f = '/data/neuralabc/paujhe/mMPI/results/4mm_flexi_ownbinmask_09062022/updatedroi_0006_D.nii.gz'
  571. tfce_sig_results = '/data/neuralabc/paujhe/mMPI/results/4mm_flexi_ownbinmask_09062022/tfce_0006_d5mored2.nii'
  573. # %%
  574. ##loop for gm_seg_f, wm_seg_f##
  575. ## follow below for d1, d2,d5 for LRN and SMP
  576. dict = {}
  577. res = {}
  578. for i in range(len(LRN)):
  579. gm_seg_f = Pname_LRN_GM_d1[i]
  580. wm_seg_f = Pname_LRN_WM_d1[i]
  581. metric_f = Pname_LRN_T1_d1[i]
  582. res[LRN[i]] = metric_extract_ROI_gm_wm_gmwm_cut(metric_f, original_ROI_f, gm_seg_f, wm_seg_f, tfce_sig_results,tfce_cut=1.3,
  583. orig_ROI_dil_iterations=3,ext_WM_dil_iterations=2,cut=0, verbosity=1)
  586. # %%
  587. df = pd.DataFrame(res)
  588. df = df.T
  589. #df.to_csv('/data/neuralabc/paujhe/mMPI/Bx/Cluster_D/august01_2024/T1mean_d1_LRN_clusterD_august01_2024.csv', index=True)
  591. # %% [markdown]
  592. # ## to check relationship between T1 values and VBM values in SPC for d1,d2,d5
  594. # %%
  595. ##### T1 values of the GM region to see associationw ith Bx
  596. dfT1_LRN_D = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/Cluster_D/august01_2024/clusterD_T1_Bx.csv') #, usecols=['T1_GM_LRN_d2','T1_GM_LRN_d5'])
  598. # %%
  599. dfT1_LRN_D['d5minusd2'] = dfT1_LRN_D['T1_GM_LRN_d5']-dfT1_LRN_D['T1_GM_LRN_d2'] #cluster D- slow learning
  600. new_dfT1_LRN_D = dfT1_LRN_D['d5minusd2']
  601. var_vec_d_T1 = new_dfT1_LRN_D.copy()
  602. var_vec_d_T1 = var_vec_d_T1[~np.isnan(var_vec_d_T1)]
  603. pearsonr(var_vec_d_T1, synB12minusB3['B12minusB3'])
  605. # %%
  606. #cluster D- left SPC
  607. T1_D_LRN_d1 = dfT1_LRN_D['T1_GM_LRN_d1']
  608. T1_D_LRN_d2 = dfT1_LRN_D['T1_GM_LRN_d2']
  609. T1_D_LRN_d5 = dfT1_LRN_D['T1_GM_LRN_d5']
  611. T1_D_SMP_d1 = dfT1_LRN_D['T1_GM_SMP_d1']
  612. T1_D_SMP_d2 = dfT1_LRN_D['T1_GM_SMP_d2']
  613. T1_D_SMP_d5 = dfT1_LRN_D['T1_GM_SMP_d5']
  615. dfvbm_D = pd.read_csv('/data/neuralabc/paujhe/mMPI/Bx/clusterD_VBM.csv')
  617. VBM_D_LRN_d1 = dfvbm_D['LRN_d1']
  618. VBM_D_LRN_d2 = dfvbm_D['LRN_d2']
  619. VBM_D_LRN_d5 = dfvbm_D['LRN_d5']
  621. VBM_D_SMP_d1 = dfvbm_D['SMP_d1']
  622. VBM_D_SMP_d2 = dfvbm_D['SMP_d2']
  623. VBM_D_SMP_d5 = dfvbm_D['SMP_d5']
  625. days = ['d5-d1','d5-d2']
  626. modality = ['gm','t1']
  627. groups = ['lrn','smp']
  629. gm_lrn_D_d1 = VBM_D_LRN_d1.copy()
  630. gm_lrn_D_d1 = gm_lrn_D_d1[~np.isnan(gm_lrn_D_d1)]
  631. gm_lrn_D_d2 = VBM_D_LRN_d2.copy()
  632. gm_lrn_D_d2 = gm_lrn_D_d2[~np.isnan(gm_lrn_D_d2)]
  633. gm_lrn_D_d5 = VBM_D_LRN_d5.copy()
  634. gm_lrn_D_d5 = gm_lrn_D_d5[~np.isnan(gm_lrn_D_d5)]
  636. gm_smp_D_d1 = VBM_D_SMP_d1.copy()
  637. gm_smp_D_d1 = gm_smp_D_d1[~np.isnan(gm_smp_D_d1)]
  638. gm_smp_D_d2 = VBM_D_SMP_d2.copy()
  639. gm_smp_D_d2 = gm_smp_D_d2[~np.isnan(gm_smp_D_d2)]
  640. gm_smp_D_d5 = VBM_D_SMP_d5.copy()
  641. gm_smp_D_d5 = gm_smp_D_d5[~np.isnan(gm_smp_D_d5)]
  643. t1_lrn_D_d1 = T1_D_LRN_d1.copy()
  644. t1_lrn_D_d1 = t1_lrn_D_d1[~np.isnan(t1_lrn_D_d1)]
  645. t1_lrn_D_d2 = T1_D_LRN_d2.copy()
  646. t1_lrn_D_d2 = t1_lrn_D_d2[~np.isnan(t1_lrn_D_d2)]
  647. t1_lrn_D_d5 = T1_D_LRN_d5.copy()
  648. t1_lrn_D_d5 = t1_lrn_D_d5[~np.isnan(t1_lrn_D_d5)]
  650. t1_smp_D_d1 = T1_D_SMP_d1.copy()
  651. t1_smp_D_d1 = t1_smp_D_d1[~np.isnan(t1_smp_D_d1)]
  652. t1_smp_D_d2 = T1_D_SMP_d2.copy()
  653. t1_smp_D_d2 = t1_smp_D_d2[~np.isnan(t1_smp_D_d2)]
  654. t1_smp_D_d5 = T1_D_SMP_d5.copy()
  655. t1_smp_D_d5 = t1_smp_D_d5[~np.isnan(t1_smp_D_d5)]
  657. gm_lrn_change = gm_lrn_D_d5 - gm_lrn_D_d2
  658. t1_lrn_change = t1_lrn_D_d5 - t1_lrn_D_d2
  659. min_len = min(len(gm_lrn_change), len(t1_lrn_change))
  660. gm_lrn_change = gm_lrn_change[:min_len]
  661. t1_lrn_change = t1_lrn_change[:min_len]
  662. r_lrn, p_lrn = pearsonr(gm_lrn_change, t1_lrn_change)
  663. ci_lower_lrn, ci_upper_lrn = correlation_ci(r_lrn, min_len)
  664. print(f"\n LRN Group (d5-d2 slow learning):")
  665. print(f" r = {r_lrn:.2f}, 95% CI [{ci_lower_lrn:.2f}, {ci_upper_lrn:.2f}], p = {p_lrn:.4f}, n = {min_len}")
  667. gm_smp_change = gm_smp_D_d5 - gm_smp_D_d2
  668. t1_smp_change = t1_smp_D_d5 - t1_smp_D_d2
  669. min_len_smp = min(len(gm_smp_change), len(t1_smp_change))
  670. gm_smp_change = gm_smp_change[:min_len_smp]
  671. t1_smp_change = t1_smp_change[:min_len_smp]
  672. r_smp, p_smp = pearsonr(gm_smp_change, t1_smp_change)
  673. ci_lower_smp, ci_upper_smp = correlation_ci(r_smp, min_len_smp)
  674. print(f"\nSMP Group (d5-d2):")
  675. print(f" r = {r_smp:.2f}, 95% CI [{ci_lower_smp:.2f}, {ci_upper_smp:.2f}], p = {p_smp:.4f}, n = {min_len_smp}")
  677. # T1 change vs SYN change for LRN group (slow learning, d5-d2)
  678. bx_data = synB12minusB3['B12minusB3'].dropna().values
  679. t1_change_bx = t1_lrn_D_d5 - t1_lrn_D_d2
  680. t1_change_bx = t1_change_bx[~np.isnan(t1_change_bx)]
  681. min_len_bx = min(len(t1_change_bx), len(bx_data))
  682. t1_change_bx = t1_change_bx[:min_len_bx]
  683. bx_data = bx_data[:min_len_bx]
  684. r_bx, p_bx = pearsonr(t1_change_bx, bx_data)
  685. ci_lower_bx, ci_upper_bx = correlation_ci(r_bx, min_len_bx)
  686. print(f"\nLeft SPC - T1 change vs SYN change (LRN, slow learning d5-d2):")
  687. print(f" r = {r_bx:.2f}, 95% CI [{ci_lower_bx:.2f}, {ci_upper_bx:.2f}], p = {p_bx:.4f}, n = {min_len_bx}")
  689. # T1 change vs SYN change for SMP group (slow learning, d5-d2) in Left SPC
  690. SMP_bx_data = synscores_SMP['B12minusB3'].dropna().values
  691. SMP_t1_change_bx = t1_smp_D_d5 - t1_smp_D_d2
  692. SMP_t1_change_bx = SMP_t1_change_bx[~np.isnan(SMP_t1_change_bx)]
  693. SMP_min_len_bx = min(len(SMP_t1_change_bx), len(SMP_bx_data))
  694. SMP_t1_change_bx = SMP_t1_change_bx[:SMP_min_len_bx]
  695. SMP_bx_data = SMP_bx_data[:SMP_min_len_bx]
  696. SMP_r_bx, SMP_p_bx = pearsonr(SMP_t1_change_bx, SMP_bx_data)
  697. SMP_ci_lower_bx, SMP_ci_upper_bx = correlation_ci(SMP_r_bx, SMP_min_len_bx)
  698. print(f"\nLeft SPC - T1 change vs SYN change (SMP, slow learning d5-d2):")
  699. print(f" r = {SMP_r_bx:.2f}, 95% CI [{SMP_ci_lower_bx:.2f}, {SMP_ci_upper_bx:.2f}], p = {SMP_p_bx:.4f}, n = {SMP_min_len_bx}")
  701. # %% [markdown]
  702. # # LOSO (LEAVE ONE SUBJECT OUT )
  704. # %%
  705. print("DIAGNOSTIC: Checking data dimensions")
  706. print(f"VBM file length: {len(dfvbm_LRN_D)}")
  707. print(f"T1 file length: {len(dfT1_LRN_D)}")
  708. print(f"Behavior file length: {len(synscores_LRN)}")
  709. print("\nChecking for missing values in each file:")
  710. print(f"VBM NaNs:\n{dfvbm_LRN_D.isnull().sum()}")
  711. print(f"\nT1 NaNs:\n{dfT1_LRN_D.isnull().sum()}")
  712. print(f"\nBehavior NaNs:\n{synscores_LRN.isnull().sum()}")
  713. if not (len(dfvbm_LRN_D) == len(dfT1_LRN_D) == len(synscores_LRN)):
  714. print("\n WARNING: Files have different lengths!")
  715. print("This will cause misalignment. Please check your CSV files.")
  716. min_len = min(len(dfvbm_LRN_D), len(dfT1_LRN_D), len(synscores_LRN))
  717. print(f"\nTruncating all to shortest length: {min_len}")
  718. dfvbm_LRN_D = dfvbm_LRN_D.iloc[:min_len]
  719. #dfT1_LRN_D = dfT1_LRN_D.iloc[:min_len]
  720. synscores_LRN = synscores_LRN.iloc[:min_len]
  722. data = pd.DataFrame({
  723. # VBM GMV values
  724. 'gmv_d1': dfvbm_LRN_D['LRN_d1'],
  725. 'gmv_d2': dfvbm_LRN_D['LRN_d2'],
  726. 'gmv_d5': dfvbm_LRN_D['LRN_d5'],
  727. 'gmv_d6': dfvbm_LRN_D['LRN_d6'],
  728. 'gmv_d1_SMP': dfvbm_SMP_D['SMP_d1'],
  729. 'gmv_d2_SMP': dfvbm_SMP_D['SMP_d2'],
  730. 'gmv_d5_SMP': dfvbm_SMP_D['SMP_d5'],
  731. 'gmv_d6_SMP': dfvbm_SMP_D['SMP_d6'],
  732. 't1_gmv_d1': dfT1_LRN_D['T1_GM_LRN_d1'],
  733. 't1_gmv_d2': dfT1_LRN_D['T1_GM_LRN_d2'],
  734. 't1_gmv_d5': dfT1_LRN_D['T1_GM_LRN_d5'],
  735. 't1_gmv_d1_SMP': dfT1_LRN_D['T1_GM_SMP_d1'],
  736. 't1_gmv_d2_SMP': dfT1_LRN_D['T1_GM_SMP_d2'],
  737. 't1_gmv_d5_SMP': dfT1_LRN_D['T1_GM_SMP_d5'],
  738. 'syn_b1': synscores_LRN['Block1'],
  739. 'syn_b3': synscores_LRN['Block3'],
  740. 'syn_b12': synscores_LRN['Block12']
  742. })
  744. # %%
  745. # SMP GROUP: T1 GMV Change (D5-D2) vs Behavior (Block12-Block3)
  746. data_SMP = pd.DataFrame({
  747. 't1_gmv_d2_SMP': dfT1_LRN_D['T1_GM_SMP_d2'],
  748. 't1_gmv_d5_SMP': dfT1_LRN_D['T1_GM_SMP_d5'],
  749. 'SMP_syn_b3': synscores_SMP['Block 3'],
  750. 'SMP_syn_b12': synscores_SMP['Block 12']
  751. }).reset_index(drop=True)
  753. print(f"SMP group n before dropna: {len(data_SMP)}")
  754. print(data_SMP.isnull().sum())
  755. data_SMP = data_SMP.dropna()
  756. print(f"SMP group n after dropna: {len(data_SMP)}")
  758. # Change scores
  759. data_SMP['t1_change_SMP'] = data_SMP['t1_gmv_d5_SMP'] - data_SMP['t1_gmv_d2_SMP']
  760. data_SMP['bx_change_SMP'] = data_SMP['SMP_syn_b12'] - data_SMP['SMP_syn_b3']
  762. X_t1_SMP = data_SMP['t1_change_SMP'].values
  763. y_SMP = data_SMP['bx_change_SMP'].values
  765. # --- Biased (standard) Pearson ---
  766. r_SMP, p_SMP = stats.pearsonr(X_t1_SMP, y_SMP)
  767. print(f"\nSMP: T1 GMV Change (D5-D2) vs Behavior Change (B12-B3)")
  768. print(f" Original (BIASED): r = {r_SMP:.3f}, p = {p_SMP:.4f}, n = {len(X_t1_SMP)}")
  770. # --- LOSO (unbiased) ---
  771. loo = LeaveOneOut()
  772. loso_SMP = []
  773. for train_idx, test_idx in loo.split(X_t1_SMP):
  774. r_train, _ = stats.pearsonr(X_t1_SMP[train_idx], y_SMP[train_idx])
  775. loso_SMP.append(r_train)
  777. mean_r_SMP = np.mean(loso_SMP)
  778. ci_lower_SMP = np.percentile(loso_SMP, 2.5)
  779. ci_upper_SMP = np.percentile(loso_SMP, 97.5)
  781. print(f" LOSO (UNBIASED): mean r = {mean_r_SMP:.3f}")
  782. print(f" 95% CI: [{ci_lower_SMP:.3f}, {ci_upper_SMP:.3f}]")
  783. print(f" SD = {np.std(loso_SMP):.3f}")
  785. # %%
  786. # Calculate changes (slow learning stage)
  787. data['gmv_change'] = data['gmv_d5'] - data['gmv_d2']
  788. data['t1_change'] = data['t1_gmv_d5'] - data['t1_gmv_d2']
  789. data['behavior_change'] = data['syn_b12'] - data['syn_b3']
  790. data['t1_change_SMP'] = data['t1_gmv_d5_SMP'] - data['t1_gmv_d2_SMP']
  791. data['gmv_change_SMP'] = data['gmv_d5_SMP'] - data['gmv_d2_SMP']
  793. if data.isnull().any().any():
  794. print("\n WARNING: Missing values detected")
  795. print(data.isnull().sum())
  796. data = data.dropna()
  797. print(f"After removing missing: n = {len(data)}")
  799. print("CORRELATION 1: VBM GMV Change (D5-D2) vs Behavior (Block12-Block3)")
  800. X_gmv_bx = data['gmv_change'].values
  801. y_bx = data['behavior_change'].values
  802. r_gmv_bx, p_gmv_bx = stats.pearsonr(X_gmv_bx, y_bx)
  803. print(f"\n Original (BIASED): r = {r_gmv_bx:.3f}, p = {p_gmv_bx:.4f}")
  805. # LOSO
  806. loo = LeaveOneOut()
  807. loso_gmv_bx = []
  808. for train_idx, test_idx in loo.split(X_gmv_bx):
  809. r_train, _ = stats.pearsonr(X_gmv_bx[train_idx], y_bx[train_idx])
  810. loso_gmv_bx.append(r_train)
  811. mean_r_gmv_bx = np.mean(loso_gmv_bx)
  812. ci_lower_gmv_bx = np.percentile(loso_gmv_bx, 2.5)
  813. ci_upper_gmv_bx = np.percentile(loso_gmv_bx, 97.5)
  815. print(f"\n LOSO (UNBIASED): mean r = {mean_r_gmv_bx:.3f}")
  816. print(f" 95% CI: [{ci_lower_gmv_bx:.3f}, {ci_upper_gmv_bx:.3f}]")
  817. print(f" SD = {np.std(loso_gmv_bx):.3f}")
  820. print("CORRELATION 2: T1 GMV Change vs VBM GMV Change (both D5-D2)")
  821. X_t1 = data['t1_change'].values
  822. X_gmv = data['gmv_change'].values
  823. r_t1_gmv, p_t1_gmv = stats.pearsonr(X_t1, X_gmv)
  824. print(f"\n Original (BIASED): r = {r_t1_gmv:.3f}, p = {p_t1_gmv:.4f}")
  826. # LOSO
  827. loso_t1_gmv = []
  828. for train_idx, test_idx in loo.split(X_t1):
  829. r_train, _ = stats.pearsonr(X_t1[train_idx], X_gmv[train_idx])
  830. loso_t1_gmv.append(r_train)
  831. mean_r_t1_gmv = np.mean(loso_t1_gmv)
  832. ci_lower_t1_gmv = np.percentile(loso_t1_gmv, 2.5)
  833. ci_upper_t1_gmv = np.percentile(loso_t1_gmv, 97.5)
  835. print(f"\n LOSO (UNBIASED): mean r = {mean_r_t1_gmv:.3f}")
  836. print(f" 95% CI: [{ci_lower_t1_gmv:.3f}, {ci_upper_t1_gmv:.3f}]")
  837. print(f" SD = {np.std(loso_t1_gmv):.3f}")
  840. print("SMP T1 GMV Change vs VBM GMV Change (both D5-D2)")
  841. X_t1_SMP = data['t1_change_SMP'].values
  842. X_gmv_SMP = data['gmv_change_SMP'].values
  843. r_t1_gmv_SMP, p_t1_gmv_SMP = stats.pearsonr(X_t1_SMP, X_gmv_SMP)
  844. print(f"\n Original (BIASED): r = {r_t1_gmv_SMP:.3f}, p = {p_t1_gmv_SMP:.4f}")
  845. # LOSO
  846. loso_t1_gmv_SMP = []
  847. for train_idx, test_idx in loo.split(X_t1_SMP):
  848. r_train_SMP, _ = stats.pearsonr(X_t1_SMP[train_idx], X_gmv_SMP[train_idx])
  849. loso_t1_gmv_SMP.append(r_train_SMP)
  850. mean_r_t1_gmv_SMP = np.mean(loso_t1_gmv_SMP)
  851. ci_lower_t1_gmv_SMP = np.percentile(loso_t1_gmv_SMP, 2.5)
  852. ci_upper_t1_gmv_SMP = np.percentile(loso_t1_gmv_SMP, 97.5)
  853. print(f"\n LOSO (UNBIASED): mean r = {mean_r_t1_gmv_SMP:.3f}")
  854. print(f" 95% CI: [{ci_lower_t1_gmv_SMP:.3f}, {ci_upper_t1_gmv_SMP:.3f}]")
  855. print(f" SD = {np.std(loso_t1_gmv_SMP):.3f}")
  857. print("CORRELATION 3: T1 GMV Change (D5-D2) vs Behavior (Block12-Block3)")
  858. X_t1_bx = data['t1_change'].values
  859. y_bx_t1 = data['behavior_change'].values
  860. r_t1_bx, p_t1_bx = stats.pearsonr(X_t1_bx, y_bx_t1)
  861. print(f"\n Original (BIASED): r = {r_t1_bx:.3f}, p = {p_t1_bx:.4f}")
  862. # LOSO
  863. loso_t1_bx = []
  864. for train_idx, test_idx in loo.split(X_t1_bx):
  865. r_train, _ = stats.pearsonr(X_t1_bx[train_idx], y_bx_t1[train_idx])
  866. loso_t1_bx.append(r_train)
  868. mean_r_t1_bx = np.mean(loso_t1_bx)
  869. ci_lower_t1_bx = np.percentile(loso_t1_bx, 2.5)
  870. ci_upper_t1_bx = np.percentile(loso_t1_bx, 97.5)
  872. print(f"\n LOSO (UNBIASED): mean r = {mean_r_t1_bx:.3f}")
  873. print(f" 95% CI: [{ci_lower_t1_bx:.3f}, {ci_upper_t1_bx:.3f}]")
  874. print(f" SD = {np.std(loso_t1_bx):.3f}")

ROI_Bx_T1_correlation.ipynb at commit dc8b802, no license · at the source

Overview

Authors: Jhelum Paul1,2, A. T. P. Jäger3,4,5, J. Huck6, C. L. Tardif7,8,9, A. Villringer5,10, C. J. Gauthier2,11,12, P. L. Bazin13, Christopher J. Steele1,2,5
ORCID iDs: Jhelum Paul, A. T. P. Jäger, C. L. Tardif, A. Villringer, C. J. Gauthier, Christopher J. Steele
13 affiliations
  1. Department of Psychology Concordia University Montreal Québec Canada
  2. School of Health Concordia University Montréal Québec Canada
  3. Brain Language Lab Freie Universität Berlin Berlin Germany
  4. Charité Universitätsmedizin Berlin Germany
  5. Department of Neurology Max Planck Institute for Human Cognitive and Brain Sciences Leipzig Germany
  6. Department of Nuclear Medicine and Radiobiology Université de Sherbrooke Sherbrooke Québec Canada
  7. Department of Biomedical Engineering McGill University Montreal Québec Canada
  8. McConnell Brain Imaging Centre Montreal Neurological Institute Montreal Québec Canada
  9. Department of Neurology and Neurosurgery McGill University Montreal Québec Canada
  10. Clinic for Cognitive Neurology Leipzig Germany
  11. Department of Physics Concordia University Montreal Québec Canada
  12. Montreal Heart Institute Montreal Québec Canada
  13. Full Brain Picture Analytics Leidein the Netherlands
Institutions: Concordia University (Canada); Max Planck Institute for Human Cognitive and Brain Sciences (Germany); Freie Universität Berlin (Germany); Charité - Universitätsmedizin Berlin (Germany); Université de Sherbrooke (Canada); Montreal Neurological Institute and Hospital (Canada); McGill University (Canada); Montreal Heart Institute (Canada)
Journal: Human brain mapping, volume 47, issue 8, article e70562
Dates: received 12 July 2025; accepted 21 May 2026; published online 4 June 2026; in print June 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1002/hbm.70562 · PMID 42237744 · PMCID PMC13266411 · OpenAlex W7163567069
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism)
Methods: Statistics, Machine learning, Preprocessing, Connectivity, fMRI & imaging
MeSH: Brain*, Gray Matter*, Learning*, Neuronal Plasticity*, Serial Learning*, Adult, Brain Mapping, Female, Humans, Image Processing, Computer-Assisted, Magnetic Resonance Imaging, Male, Young Adult (* major topic)
Topic: Transcranial Magnetic Stimulation Studies (Neurology, Neuroscience), according to OpenAlex
Funding: Natural Sciences and Engineering Research Council of Canada (RGPIN‐2020‐06812, DGECR‐2020‐00146); Heart and Stroke Foundation of Canada New Investigator Award and Catalyst Award from the Canadian Institutes of Health Research (HNC 170723); Fonds de Recherche du Québec Santé Chercheurs‐boursier (FRQS CB Junior 2 349443)
Citations: not cited yet (Europe PMC); 119 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.

neuralabc/Structural_Physiology_MSLplasticity

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: dc8b8028c55c9a2d24fb99c44173d520de00cf1e, 15 April 2026
Languages: MATLAB (5), Jupyter (1)
Size: 10 files, 6 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: SPM (5 files), Matplotlib (1 file), MRtrix3 (1 file), NiBabel (1 file), Nilearn (1 file), NumPy (1 file), pandas (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
7 files
At the source: github.com/neuralabc/Structural_Physiology_MSLplasticity

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 6 scripts, each with its path and the digest of its content;
  • 7 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1002/hbm.70562.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 13 MeSH terms, 3 funders, 117 references.

Cite

This paper

Paul, J., Jäger, A. T. P., Huck, J., Tardif, C. L., Villringer, A., Gauthier, C. J., Bazin, P. L., & Steele, C. J. (2026). Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning. Human brain mapping, 47(8), e70562. https://doi.org/10.1002/hbm.70562

BibTeX

@article{paul2026distinct,
author = {Paul, Jhelum and Jäger, A. T. P. and Huck, J. and Tardif, C. L. and Villringer, A. and Gauthier, C. J. and Bazin, P. L. and Steele, Christopher J.},
title = {{Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning}},
journal = {Human brain mapping},
year = {2026},
month = jun,
volume = {47},
number = {8},
pages = {e70562},
publisher = {Wiley},
issn = {1065-9471},
doi = {10.1002/hbm.70562},
url = {https://doi.org/10.1002/hbm.70562},
pmid = {42237744},
pmcid = {PMC13266411}
}

RIS

TY - JOUR
AU - Paul, Jhelum
AU - Jäger, A. T. P.
AU - Huck, J.
AU - Tardif, C. L.
AU - Villringer, A.
AU - Gauthier, C. J.
AU - Bazin, P. L.
AU - Steele, Christopher J.
TI - Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning
T2 - Human brain mapping
J2 - Hum Brain Mapp
PY - 2026
DA - 2026/06/01
VL - 47
IS - 8
SP - e70562
SN - 1065-9471
PB - Wiley
DO - 10.1002/hbm.70562
UR - https://doi.org/10.1002/hbm.70562
LA - en
ER -

CSL-JSON

{
"id": "10.1002/hbm.70562",
"type": "article-journal",
"title": "Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning",
"container-title": "Human brain mapping",
"author": [
{
"family": "Paul",
"given": "Jhelum"
},
{
"family": "Jäger",
"given": "A. T. P."
},
{
"family": "Huck",
"given": "J."
},
{
"family": "Tardif",
"given": "C. L."
},
{
"family": "Villringer",
"given": "A."
},
{
"family": "Gauthier",
"given": "C. J."
},
{
"family": "Bazin",
"given": "P. L."
},
{
"family": "Steele",
"given": "Christopher J."
}
],
"container-title-short": "Hum Brain Mapp",
"volume": "47",
"issue": "8",
"page": "e70562",
"DOI": "10.1002/hbm.70562",
"PMID": "42237744",
"PMCID": "PMC13266411",
"ISSN": "1065-9471",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/hbm.70562",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-73366-9 [code]
Cortical and white matter myelination proceed in concert during early infancy.
Journal: Nature communications
In common: MRtrix3, statsmodels, NiBabel, 6 other tools, 6 references
[2] doi:10.1162/imag.a.1325 [code]
Decoding everyday levels of musical training from subcortical white-matter architecture.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: MRtrix3, statsmodels, seaborn, 4 other tools, structural MRI / diffusion, 6 references
[3] doi:10.1162/imag.a.1203 [code]
Motor cortical areas facilitate schema-mediated integration of new motor information into memory.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: SPM, 8 references
[4] doi:10.1162/imag.a.1174 [code]
Retinal nerve fibre layer thickness reflects characteristics of brain grey and white matter.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, SPM, NiBabel, 2 other tools, structural MRI / diffusion, 2 references, author A. Villringer
[5] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: MRtrix3, Nilearn, SPM, 8 other tools
[6] doi:10.1038/s41467-026-71918-7 [code]
Developmental disinhibition gates language lateralization in childhood.
Journal: Nature communications
In common: MRtrix3, Nilearn, SPM, 7 other tools, 1 reference
[7] doi:10.1162/imag.a.105 [code]
Right posterior theta reflects human parahippocampal phase resetting by salient cues during goal-directed navigation
Journal: n/a
In common: Nilearn, SPM, statsmodels, 6 other tools, 2 references
[8] doi:10.1038/s41467-026-74566-z [code]
Low-dimensional and optimised representations of high-level information in the expert brain.
Journal: Nature communications
In common: Nilearn, SPM, statsmodels, 7 other tools, 1 reference
[9] doi:10.1038/s41467-026-71568-9 [code]
Convergent and selective representations of pain, appetitive processes, aversive processes, and cognitive control in the insula.
Journal: Nature communications
In common: MRtrix3, Nilearn, SPM, 7 other tools
[10] doi:10.7554/elife.107933 [code]
Modality-agnostic decoding of vision and language from fMRI.
Journal: eLife
In common: Nilearn, SPM, statsmodels, 7 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.