OSCR

Cortical and white matter myelination proceed in concert during early infancy.

Code ↔ Paper

11 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 11 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Bundle identification with pyBabyAFQ ↔ pyafq_pipeline/pyafq_pipeline_dHCP.ipynb, the whole file · a weak match · score 0.74 · babyAFQ, waypoints, template, AWS, MRtrix, space
  2. [2] § Methods › Assessing myelin-sensitive imaging metrics across tissues ↔ mrtrix_pipeline/mrtrix_pipeline_dHCP.ipynb, lines 345–389 · score 0.73 · draw EM, FreeSurfer, definition, boundaries, MRtrix, segmentations
  3. [3] § Methods › Tractography ↔ mrtrix_pipeline/mrtrix_pipeline_dHCP.ipynb, lines 442–504 · score 0.65 · IFOD1, Seeds, angle, dHCP, MRtrix, CSD
  4. [4] § Methods › Statistical analyses ↔ Figure4.ipynb, lines 217–240 · score 0.61 · motor scores, birth age, ANOVA, scan age, sex, model
  5. [5] § Results › Coupled myelin development of bundles and their cortical targets ↔ Figure3.ipynb, lines 409–543 · score 0.59 · FcMa, FcMi, regression, hues, IFOF, MLF
  6. [6] § Results › Reduced T1w/T2w coupling in preterm infants ↔ Figure5.ipynb, lines 154–192 · score 0.55 · FcMa, FcMi, IFOF, MLF, AF, UNC
  7. [7] § Methods › Tractography ↔ pyafq_pipeline/pyafq_pipeline_dHCP.ipynb, the whole file · a weak match · score 0.55 · tracking, MRtrix, CSD, tractography, segmentation, streamlines
  8. [8] § Results › Coupled myelin content of bundles and their cortical targets ↔ Figure3.ipynb, lines 409–543 · score 0.54 · FcMa, FcMi, regression, IFOF, MLF, AF
  9. [9] § Methods › Bundle identification with pyBabyAFQ ↔ EndpointsToCortex.ipynb, lines 632–719 · score 0.53 · Python, ROIs, frontal, cingulate, occipital, inferior
  10. [10] § Results › Coupled myelin content of bundles and their cortical targets ↔ Figure2.ipynb, lines 21–103 · score 0.51 · FcMa, FcMi, IFOF, MLF, AF, UNC
  11. [11] § Methods › Statistical analyses ↔ Figure3.ipynb, lines 213–309 · score 0.50 · birth age, scan age, variables, predicting, model, linear

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 · 1,128 lines · 39 KB · no license · 3 matches

  1. # %% [markdown]
  2. # ## Set up environment
  4. # %%
  5. # Set up environment
  6. import warnings;
  7. warnings.filterwarnings('ignore');
  8. from IPython.core.interactiveshell import InteractiveShell;
  9. InteractiveShell.ast_node_interactivity = "all";
  10. import pandas as pd
  11. import os
  12. import numpy as np
  13. import seaborn as sns
  14. import matplotlib.pyplot as plt
  15. import scipy
  16. from scipy import stats
  17. from matplotlib.lines import Line2D;
  18. from sklearn import preprocessing
  19. import statsmodels.formula.api as smf
  20. from utils import get_data
  22. # %%
  23. #plotting parameters, these will be the same for all plots
  25. bundles = ['AFL', 'AFR', 'ATRL', 'ATRR', 'CCL', 'CCR', 'CSL', 'CSR', 'FcMi', 'FcMa', 'IFOFL', 'IFOFR', 'ILFL', 'ILFR',
  26. 'MLFL', 'MLFR', 'ORL', 'ORR', 'SLFL', 'SLFR', 'UNCL', 'UNCR', 'VOFL', 'VOFR', 'pAFL', 'pAFR']
  27. rightBundles = ['AFR', 'ATRR', 'CCR', 'CSR', 'FcMi', 'FcMa', 'IFOFR', 'ILFR', 'MLFR', 'ORR', 'pAFR', 'SLFR', 'UNCR', 'VOFR']
  28. leftBundles = ['AFL', 'ATRL', 'CCL', 'CSL', 'FcMi', 'FcMa', 'IFOFL', 'ILFL', 'MLFL', 'ORL', 'pAFL', 'SLFL', 'UNCL', 'VOFL']
  31. tractPosHorz = {'AF': (0,0), 'ATR': (0, 1), 'CC': (0, 2), 'CS': (0, 3),
  32. 'IFOF':(1,0), 'ILF':(1,1), 'MLF':(1,2), 'OR': (1,3),
  33. 'SLF': (2, 0), 'UNC':(2,1), 'VOF': (2, 2), 'pAF': (2,3),
  34. 'FcMi': (3, 1), 'FcMa': (3, 2)};
  36. tracts=['AF', 'ATR', 'CC', 'CS', 'IFOF', 'ILF', 'MLF', 'OR', 'SLF', 'UNC', 'VOF', 'pAF', 'FcMi', 'FcMa']
  37. colors=['cyan', 'blue', 'green', 'orange', 'purple', 'brown', 'olive', 'coral', 'fuchsia', 'yellow', 'indigo', 'violet', 'salmon', 'red']
  39. color_list_all=sns.color_palette("tab20")+sns.color_palette("tab20b")
  40. color_list_complete=color_list_all
  41. color_list_chosen=color_list_all[18:20]
  42. color_list_chosen.extend(color_list_all[0:2])
  43. color_list_chosen.extend(color_list_all[4:6])
  44. color_list_chosen.extend(color_list_all[2:4])
  45. color_list_chosen.extend(color_list_all[7:8])
  46. color_list_chosen.extend(color_list_all[6:7])
  47. color_list_chosen.extend(color_list_all[8:10])
  48. color_list_chosen.extend(color_list_all[10:12])
  49. color_list_chosen.extend(color_list_all[24:25])
  50. color_list_chosen.extend(color_list_all[27:28])
  51. color_list_chosen.extend(color_list_all[34:35])
  52. color_list_chosen.extend(color_list_all[35:36])
  53. color_list_chosen.extend(color_list_all[12:14])
  54. color_list_chosen.extend(color_list_all[16:18])
  55. color_list_chosen.extend(color_list_all[22:23])
  56. color_list_chosen.extend(color_list_all[23:24])
  57. color_list_chosen.extend(color_list_all[37:38])
  58. color_list_chosen.extend(color_list_all[39:40])
  60. sns.palplot(color_list_all)
  61. sns.palplot(color_list_chosen)
  62. color_list_all=color_list_chosen
  63. color_order=[19, 18, 1, 0, 5, 4, 3, 2, 7, 6, 9, 8, 11, 10, 15, 14, 21, 20, 25, 24, 13, 12, 17, 16, 23, 22]
  65. color_list_nohemis=sns.color_palette("tab20")+sns.color_palette("tab20b")
  66. color_list_complete_nohemis=color_list_nohemis
  67. color_list_chosen_nohemis=color_list_nohemis[18:19]
  68. color_list_chosen_nohemis.extend(color_list_nohemis[0:1])
  69. color_list_chosen_nohemis.extend(color_list_nohemis[4:5])
  70. color_list_chosen_nohemis.extend(color_list_nohemis[2:3])
  71. color_list_chosen_nohemis.extend(color_list_nohemis[8:9])
  72. color_list_chosen_nohemis.extend(color_list_nohemis[10:11])
  73. color_list_chosen_nohemis.extend(color_list_nohemis[24:25])
  74. color_list_chosen_nohemis.extend(color_list_nohemis[34:35])
  75. color_list_chosen_nohemis.extend(color_list_nohemis[12:13])
  76. color_list_chosen_nohemis.extend(color_list_nohemis[16:17])
  77. color_list_chosen_nohemis.extend(color_list_nohemis[22:23])
  78. color_list_chosen_nohemis.extend(color_list_nohemis[37:38])
  79. color_list_chosen_nohemis.extend(color_list_nohemis[7:8])
  80. color_list_chosen_nohemis.extend(color_list_nohemis[6:7])
  81. sns.palplot(color_list_chosen_nohemis)
  83. color_order2=[19, 18, 5, 4, 9, 8, 11, 10, 15, 14, 21, 20, 25, 24, 13, 12, 17, 16, 23, 22, 7, 6]
  85. color_list_3=sns.color_palette("tab20")+sns.color_palette("tab20b")
  86. color_list_complete3=color_list_3
  87. color_list_chosen3=color_list_3[18:19]
  88. color_list_chosen3.extend(color_list_3[0:1])
  89. color_list_chosen3.extend(color_list_3[4:5])
  90. color_list_chosen3.extend(color_list_3[2:3])
  91. color_list_chosen3.extend(color_list_3[8:9])
  92. color_list_chosen3.extend(color_list_3[10:11])
  93. color_list_chosen3.extend(color_list_3[24:25])
  94. color_list_chosen3.extend(color_list_3[34:35])
  95. color_list_chosen3.extend(color_list_3[12:13])
  96. color_list_chosen3.extend(color_list_3[16:17])
  97. color_list_chosen3.extend(color_list_3[22:23])
  98. color_list_chosen3.extend(color_list_3[37:38])
  99. color_list_chosen3.extend(color_list_3[7:8])
  100. color_list_chosen3.extend(color_list_3[6:7])
  101. sns.palplot(color_list_chosen3)
  104. color_order3=[19, 1, 5, 3, 9, 11, 25, 35, 13, 17, 23, 38, 8, 7]
  106. # %%
  107. get_data()
  109. # %%
  110. slope_analyses_mean=pd.read_csv('./inputData/GMandWMT1wT2wSubjectsAge.csv', index_col=None)
  111. SlopeDataframe=pd.read_csv('./inputData/R1Slope.csv', index_col=None)
  113. # %% [markdown]
  114. # ## T1w/T2w
  116. # %%
  117. #Compute slope for scan age for T1w/T2w of white matter
  119. tractCount=slope_analyses_mean['tractID'].unique()
  121. coeff=[] #Creating arrays
  122. se = []
  123. coeffType=[]
  124. tractTrack1=[]
  125. intercept=[]
  126. r2 = []
  127. adjustR2 = []
  128. aic = []
  129. pvals=[]
  131. for tract in tractCount:
  133. dfForStats=slope_analyses_mean[(slope_analyses_mean['tractID']==tract)]
  134. dfForStats.reset_index(level=0, inplace=True)
  136. scanAgeArray=np.array(dfForStats['scan_age'])
  137. scanAgeArray=scanAgeArray.reshape(1,-1)
  138. scanAgeNorm=preprocessing.normalize(scanAgeArray,axis=1)
  139. scanAgeNorm=scanAgeNorm.reshape(-1)
  140. dfForStats["scanAgeNorm"]=scanAgeNorm
  142. birthAgeArray=np.array(dfForStats['birth_age'])
  143. birthAgeArray=birthAgeArray.reshape(1,-1)
  144. birthAgeNorm=preprocessing.normalize(birthAgeArray,axis=1)
  145. birthAgeNorm=birthAgeNorm.reshape(-1)
  146. dfForStats["birthAgeNorm"]=birthAgeNorm
  148. t1wt2wArray=np.array(dfForStats['t1wt2w'])
  149. t1wt2wArray=t1wt2wArray.reshape(1,-1)
  150. t1wt2wNorm=preprocessing.normalize(t1wt2wArray,axis=1)
  151. t1wt2wNorm=t1wt2wNorm.reshape(-1)
  152. dfForStats["t1wt2wNorm"]=t1wt2wNorm
  154. df = dfForStats
  155. md = smf.ols(formula="t1wt2w~ 1 + scan_age", data=df) #Predicting dti_mdNorm through the other variables
  156. mdf = md.fit() # Actually commanding the linear analysis, it could be any other analysis method
  157. results=mdf.summary()
  159. # Note that tables is a list. The table at index 1 is the "core" table. Additionally, read_html puts dfs in a list, so we want index 0
  160. results_as_html = results.tables[1].as_html()
  161. LM_results=pd.read_html(results_as_html, header=0, index_col=0)[0]
  163. results_as_html_R2 = results.tables[0].as_html()
  164. LM_results_R2=pd.read_html(results_as_html_R2, index_col=0)[0]
  166. results_R2s=LM_results_R2[3]
  167. #LM_results.head()
  168. #LM_results_R2.head()
  169. #results_R2s.head()
  171. SE=LM_results['std err']
  173. inter=mdf.params['Intercept']
  174. intercept.append(inter)
  176. #SE.head()
  177. se.append(SE["scan_age"])
  178. coeff.append(mdf.params['scan_age'])
  179. coeff_type=3
  180. tractTrack1.append(tract)
  181. #nodeTrack1.append(node)
  182. coeffType.append(coeff_type) #Add another line below
  183. r2.append(results_R2s['Dep. Variable:'])
  184. adjustR2.append(results_R2s['Model:'])
  185. aic.append(results_R2s['No. Observations:'])
  186. pvals.append(results_R2s['Date:'])
  188. print(mdf.summary())
  189. df=pd.DataFrame(tractTrack1)
  190. df.columns=['tractID']
  191. df.insert(1,"coeffType",coeffType)
  192. df.insert(2,"coeff",coeff)
  193. df.insert(3,"se",se)
  194. df.insert(4,"intercept",intercept)
  195. df.insert(5,"r2",r2)
  196. df.insert(6,"adjustR2",adjustR2)
  197. df.insert(7,"aic",aic)
  198. df.insert(8,"pvals",pvals)
  200. os.makedirs('./outputData',exist_ok=True)
  201. df.to_csv('./outputData/LM_result_scanAge_noNorm_WM_AverageSubj.csv')
  202. df.head()
  203. df.min(axis=0, numeric_only=True)
  204. df.max(axis=0, numeric_only=True)
  205. df.mean(axis=0, numeric_only=True)
  206. dfScanAge=df
  207. df.head()
  209. statsPerBundle=df.groupby('tractID').mean()
  210. statsPerBundle.head()
  211. statsPerBundle.to_csv('./outputData/LM_result_scanAge_noNorm_meanPerBundle_WM_AverageSubj.csv')
  213. # %%
  214. #Compute slope for scan age for T1w/T2w gray matter
  216. tractCount=slope_analyses_mean['tractID'].unique()
  218. coeff=[] #Creating arrays
  219. se = []
  220. coeffType=[]
  221. tractTrack1=[]
  222. intercept=[]
  223. r2 = []
  224. adjustR2 = []
  225. aic = []
  226. pvals=[]
  228. for tract in tractCount:
  230. dfForStats=slope_analyses_mean[(slope_analyses_mean['tractID']==tract)]
  231. dfForStats.reset_index(level=0, inplace=True)
  233. scanAgeArray=np.array(dfForStats['scan_age'])
  234. scanAgeArray=scanAgeArray.reshape(1,-1)
  235. scanAgeNorm=preprocessing.normalize(scanAgeArray,axis=1)
  236. scanAgeNorm=scanAgeNorm.reshape(-1)
  237. dfForStats["scanAgeNorm"]=scanAgeNorm
  239. birthAgeArray=np.array(dfForStats['birth_age'])
  240. birthAgeArray=birthAgeArray.reshape(1,-1)
  241. birthAgeNorm=preprocessing.normalize(birthAgeArray,axis=1)
  242. birthAgeNorm=birthAgeNorm.reshape(-1)
  243. dfForStats["birthAgeNorm"]=birthAgeNorm
  245. t1wt2wArray=np.array(dfForStats['WeAvGMT1wT2w'])
  246. t1wt2wArray=t1wt2wArray.reshape(1,-1)
  247. t1wt2wNorm=preprocessing.normalize(t1wt2wArray,axis=1)
  248. t1wt2wNorm=t1wt2wNorm.reshape(-1)
  249. dfForStats["t1wt2wNorm"]=t1wt2wNorm
  250. #dfForStats.head()
  252. #dfForStats.to_csv('./'+tract+'_dataCleanedForLME.csv')
  253. df = dfForStats
  254. md = smf.ols(formula="WeAvGMT1wT2w~ 1 + scan_age", data=df) #Predicting dti_mdNorm through the other variables
  255. mdf = md.fit() # Actually commanding the linear analysis, it could be any other analysis method
  256. results=mdf.summary()
  258. # Note that tables is a list. The table at index 1 is the "core" table. Additionally, read_html puts dfs in a list, so we want index 0
  259. results_as_html = results.tables[1].as_html()
  260. LM_results=pd.read_html(results_as_html, header=0, index_col=0)[0]
  262. results_as_html_R2 = results.tables[0].as_html()
  263. LM_results_R2=pd.read_html(results_as_html_R2, index_col=0)[0]
  265. results_R2s=LM_results_R2[3]
  266. #LM_results.head()
  267. #LM_results_R2.head()
  268. #results_R2s.head()
  270. SE=LM_results['std err']
  272. inter=mdf.params['Intercept']
  273. intercept.append(inter)
  275. #SE.head()
  276. se.append(SE["scan_age"])
  277. coeff.append(mdf.params['scan_age'])
  278. coeff_type=3
  279. tractTrack1.append(tract)
  280. #nodeTrack1.append(node)
  281. coeffType.append(coeff_type) #Add another line below
  282. r2.append(results_R2s['Dep. Variable:'])
  283. adjustR2.append(results_R2s['Model:'])
  284. aic.append(results_R2s['No. Observations:'])
  285. pvals.append(results_R2s['Date:'])
  287. print(mdf.summary())
  288. df=pd.DataFrame(tractTrack1)
  289. df.columns=['tractID']
  290. df.insert(1,"coeffType",coeffType)
  291. df.insert(2,"coeff_GM",coeff)
  292. df.insert(3,"se",se)
  293. df.insert(4,"intercept",intercept)
  294. df.insert(5,"r2",r2)
  295. df.insert(6,"adjustR2",adjustR2)
  296. df.insert(7,"aic",aic)
  297. df.insert(8,"pvals",pvals)
  299. df.to_csv('./outputData/LM_result_scanAge_noNorm_GM_AverageSubj.csv')
  300. df.head()
  301. df.min(axis=0, numeric_only=True)
  302. df.max(axis=0, numeric_only=True)
  303. df.mean(axis=0, numeric_only=True)
  304. dfScanAge=df
  305. df.head()
  307. statsPerBundle=df.groupby('tractID').mean()
  308. statsPerBundle.head()
  309. statsPerBundle.to_csv('./outputData/LM_result_scanAge_noNorm_meanPerBundle_GM_AverageSubj.csv')
  311. # %%
  312. #Merge statistic dataframes of gray and white matter
  314. SlopeWM = pd.read_csv('./outputData/LM_result_scanAge_noNorm_WM_AverageSubj.csv', index_col=0)
  315. SlopeGM = pd.read_csv('./outputData/LM_result_scanAge_noNorm_GM_AverageSubj.csv', index_col=0)
  317. SlopeBoth = pd.merge(SlopeWM, SlopeGM[['tractID', 'coeff_GM', 'r2', 'aic']], on='tractID')
  319. SlopeBoth=SlopeBoth.set_index(['tractID']).reindex(['AFL', 'AFR', 'ATRL', 'ATRR', 'CCL', 'CCR', 'CSL', 'CSR', 'FcMi', 'FcMa', 'IFOFL', 'IFOFR', 'ILFL', 'ILFR', 'MLFL', 'MLFR', 'ORL', 'ORR', 'SLFL', 'SLFR', 'UNCL', 'UNCR', 'VOFL', 'VOFR', 'pAFL', 'pAFR']).reset_index()
  320. SlopeBoth
  321. SlopeBoth.to_csv('./outputData/SlopeT1wT2w.csv')
  323. # %%
  324. # Paired t-test between white matter and gray matter coefficients
  325. t_stat, p_val = stats.ttest_rel(SlopeBoth['coeff'], SlopeBoth['coeff_GM'])
  326. print(f"CoeffWM vs. CoeffGM t = {t_stat:.4f}, p = {p_val:.4f}")
  328. # Print mean values
  329. mean_wm = SlopeBoth['coeff'].mean()
  330. mean_gm = SlopeBoth['coeff_GM'].mean()
  331. print(f"Mean Coeff WM: {mean_wm:.4f}")
  332. print(f"Mean Coeff GM: {mean_gm:.4f}")
  334. # Calculate difference metrics
  335. SlopeBoth['diff_r2'] = SlopeBoth['r2_x'] - SlopeBoth['r2_y']
  336. SlopeBoth['diff_aic'] = SlopeBoth['aic_x'] - SlopeBoth['aic_y']
  338. # Print mean differences
  339. mean_diff_r2 = SlopeBoth['diff_r2'].mean()
  340. mean_diff_aic = SlopeBoth['diff_aic'].mean()
  341. print(f"Mean Δr²: {mean_diff_r2:.4f}")
  342. print(f"Mean ΔAIC: {mean_diff_aic:.4f}")
  344. # %%
  345. sns.set_style('white');
  346. fig1=sns.lmplot(
  347. data=SlopeBoth, x='coeff_GM', y="coeff", hue="tractID", height=10, scatter_kws={"s": 400}, fit_reg=False, legend=False, palette=color_list_chosen
  348. )
  349. fig1=sns.regplot(data=SlopeBoth, x='coeff_GM', y="coeff", scatter=False, ax=fig1.axes[0, 0], line_kws={"color": "darkgrey"})
  350. res=scipy.stats.pearsonr(SlopeBoth['coeff_GM'], SlopeBoth['coeff'])
  351. print(res)
  352. res.confidence_interval()
  353. #fig1.legend(title='Tract', fontsize='15', title_fontsize='20', bbox_to_anchor=(0.7, 0.25, 0.5, 0.5), loc='right', borderaxespad=0, frameon=FcMilse)
  354. plt.xlabel("T1w/T2w Slope in GM", fontsize=50)
  355. plt.ylabel("T1w/T2w Slope in WM", fontsize=50)
  356. plt.xticks([0.025, 0.035, 0.045], fontsize=45)
  357. plt.yticks(fontsize=45)
  359. plt.savefig('./figures/T1wT2w_SlopeGMvsSlopeWM.png', dpi=600, bbox_inches = "tight")
  361. # %%
  362. # Define bundles
  363. b=['AFL', 'AFR', 'ATRL', 'ATRR', 'CCL', 'CCR', 'CSL', 'CSR', 'FcMi', 'FcMa', 'IFOFL', 'IFOFR', 'ILFL', 'ILFR', 'MLFL', 'MLFR', 'ORL', 'ORR', 'SLFL', 'SLFR', 'UNCL', 'UNCR', 'VOFL', 'VOFR', 'pAFL', 'pAFR']
  365. results = []
  367. for bundle in b:
  368. bundle_df = slope_analyses_mean[slope_analyses_mean['tractID'] == bundle]
  370. if len(bundle_df) >= 2:
  371. # White matter
  372. r_wm, p_wm = stats.pearsonr(bundle_df['scan_age'], bundle_df['t1wt2w'])
  374. # Gray matter
  375. r_gm, p_gm = stats.pearsonr(bundle_df['scan_age'], bundle_df['WeAvGMT1wT2w'])
  378. results.append({
  379. 'Tract': bundle,
  380. 'r_WM': r_wm,
  381. 'p_WM': p_wm,
  382. 'r2_WM': r_wm**2,
  383. 'r_GM': r_gm,
  384. 'p_GM': p_gm,
  385. 'r2_GM': r_gm**2
  386. })
  388. # Convert to DataFrame
  389. df_tract_corrs = pd.DataFrame(results)
  391. # View or export
  392. df_tract_corrs
  393. df_tract_corrs.to_csv('./outputData/tract_r2_pvalues_wm_gm_bothemis.csv', index=False)
  395. r2_min_wm = df_tract_corrs['r2_WM'].min()
  396. r2_min_gm = df_tract_corrs['r2_GM'].min()
  397. r_min_wm = df_tract_corrs['r_WM'].min()
  398. r_min_gm = df_tract_corrs['r_GM'].min()
  399. p_max_wm = df_tract_corrs['p_WM'].max()
  400. p_max_gm = df_tract_corrs['p_GM'].max()
  402. print(f"Min r² WM: {r2_min_wm}")
  403. print(f"Min r² GM: {r2_min_gm}")
  404. print(f"Min r WM: {r_min_wm}")
  405. print(f"Min r GM: {r_min_gm}")
  406. print(f"Max p WM: {p_max_wm}")
  407. print(f"Max p GM: {p_max_gm}")
  409. # %%
  410. # --- PREPARE DATA -------------------------------------------------------------
  411. df = slope_analyses_mean.copy()
  413. # Extract Hemisphere from tract names (e.g. "AF_L" → "L")
  414. df["Hemisphere"] = df["tractID"].str.extract(r'(L|R)$')
  416. # Extract base bundle name (e.g. "AF_L" → "AF")
  417. df["BundleBase"] = df["tractID"].str.replace(r'(L|R)$', '', regex=True)
  419. # --- DEFINE LAYOUT AND STYLE --------------------------------------------------
  420. tractPosHorz = {
  421. 'AF': (0,0), 'ATR': (0,1), 'CC': (0,2), 'CS': (0,3),
  422. 'IFOF':(1,0), 'ILF':(1,1), 'MLF':(1,2), 'OR':(1,3),
  423. 'SLF':(2,0), 'UNC':(2,1), 'VOF':(2,2), 'pAF':(2,3),
  424. 'FcMi':(3,1), 'FcMa':(3,2)
  425. }
  427. tracts = ['AF', 'ATR', 'CC', 'CS', 'IFOF', 'ILF', 'MLF',
  428. 'OR', 'SLF', 'UNC', 'VOF', 'pAF', 'FcMi', 'FcMa']
  430. colors = ['cyan', 'blue', 'green', 'orange', 'purple', 'brown', 'olive',
  431. 'coral', 'fuchsia', 'yellow', 'indigo', 'violet', 'salmon', 'red']
  433. sns.set(font_scale=2)
  434. sns.set_style("white")
  436. fig, axes = plt.subplots(4, 4, figsize=(14,17), frameon=False)
  437. fig.subplots_adjust(wspace=0.1, hspace=0.9)
  439. import matplotlib.colors as mcolors
  441. def lighten_color(color, amount=0.5):
  442. """
  443. Lightens the given color by blending it with white.
  444. amount=0 → white, amount=1 → original color.
  445. """
  446. try:
  447. c = mcolors.cnames[color]
  448. except KeyError:
  449. c = color
  450. rgb = mcolors.to_rgb(c)
  451. return tuple(1 - (1 - x) * amount for x in rgb)
  453. # --- MAIN LOOP ---------------------------------------------------------------
  454. for ct, (bundle, color) in enumerate(zip(tracts, colors)):
  455. if bundle not in tractPosHorz:
  456. continue
  458. ax = axes[tractPosHorz[bundle]]
  460. # Filter data for the base bundle (both hemispheres)
  461. SlopeTract = df.query("BundleBase == @bundle")
  463. # Aggregate per subject, hemisphere, etc.
  464. SlopeTractGM = SlopeTract.groupby(
  465. ['subjectID', 'scan_age', 'BundleBase', 'Hemisphere'],
  466. as_index=False, sort=False
  467. )['WeAvGMT1wT2w'].mean()
  469. SlopeTractWM = SlopeTract.groupby(
  470. ['subjectID', 'scan_age', 'BundleBase', 'Hemisphere'],
  471. as_index=False, sort=False
  472. )['t1wt2w'].mean()
  474. # Make a lighter version of the current tract color for the right hemisphere
  475. light_color = lighten_color(color, amount=0.5)
  477. # Plot white matter (same hue, lighter for R)
  478. sns.scatterplot(
  479. data=SlopeTractWM, x='scan_age', y='t1wt2w', hue='Hemisphere', style='Hemisphere',
  480. markers={'L': 'o', 'R': 'o'}, s=50,
  481. palette={'L': color, 'R': light_color}, ax=ax,
  482. alpha=0.5, legend=False
  483. )
  485. # Regression lines per hemisphere
  486. for hemi, line_color in zip(['L', 'R'], [color, light_color]):
  487. sns.regplot(
  488. data=SlopeTractWM[SlopeTractWM['Hemisphere'] == hemi],
  489. x='scan_age', y='t1wt2w', scatter=False, ax=ax,
  490. line_kws={'color': line_color, 'alpha': 0.8}
  491. )
  493. # Plot gray matter (with slightly lighter colors)
  494. sns.scatterplot(
  495. data=SlopeTractGM, x='scan_age', y='WeAvGMT1wT2w', hue='Hemisphere', style='Hemisphere',
  496. markers={'L': 'o', 'R': 'o'}, s=50,
  497. palette={'L': 'darkgray', 'R': 'lightgray'}, ax=ax,
  498. alpha=0.5, legend=False
  499. )
  501. for hemi, line_color in zip(['L', 'R'], ['darkgray', 'lightgray']):
  502. sns.regplot(
  503. data=SlopeTractGM[SlopeTractGM['Hemisphere'] == hemi],
  504. x='scan_age', y='WeAvGMT1wT2w', scatter=False, ax=ax,
  505. line_kws={'color': line_color, 'alpha': 0.8}
  506. )
  508. # --- Titles & labels
  509. ax.set_title(bundle, pad=5, fontsize=25)
  510. ax.set_xlabel('Scan Age', labelpad=15, fontsize=20)
  512. if ct in [0, 4, 8, 12]: # first column
  513. _=ax.set_ylabel("T1w/T2w",fontsize=20);
  514. _=ax.set_xticks([30, 34, 38, 42]);
  515. #_=ax.set_xticklabels(fontsize=20);
  516. #_=ax.set_yticks([0.4, 0.45, 0.5]);
  517. #_=ax.set_yticklabels(fontsize=20)
  518. _=ax.spines['right'].set_visible(False);
  519. _=ax.spines['top'].set_visible(False);
  520. _=line1=Line2D([],[],color=color_list_all[ct],linestyle='dashed');
  521. _=line2=Line2D([],[],color=color_list_all[ct],linestyle='-');
  523. else:
  524. _=ax.set_ylabel('T1w/T2w', fontsize=14);
  525. _=ax.set_xticks([30, 34, 38, 42]);
  526. #_=ax.set_xticklabels(fontsize=20);
  527. #_=ax.set_yticks([0.4, 0.45, 0.5]);
  528. _=ax.yaxis.set_visible(False);
  529. _=ax.spines['right'].set_visible(False);
  530. _=ax.spines['left'].set_visible(False);
  531. _=ax.spines['top'].set_visible(False);
  533. # Clean frame
  534. ax.spines['right'].set_visible(False)
  535. ax.spines['top'].set_visible(False)
  537. # Turn off unused axes
  538. axes[3,3].axis("off")
  539. axes[3,0].axis("off")
  541. # --- SAVE & SHOW -------------------------------------------------------------
  542. fig.savefig('./figures/scanage_per_bundle_bothhemis.png', dpi=600)
  543. plt.show()
  545. # %%
  546. # --- PREPARE DATA -------------------------------------------------------------
  547. df = slope_analyses_mean.copy()
  549. # 1) Extract Hemisphere (only if tract ends with L or R)
  550. df["Hemisphere"] = df["tractID"].str.extract(r'(L|R)$')
  552. # 2) Extract BundleBase safely
  553. df["BundleBase"] = df["tractID"] # start with full name
  554. mask = df["tractID"].str.endswith(("L", "R")) # true hemisphere bundles
  555. df.loc[mask, "BundleBase"] = df.loc[mask, "tractID"].str.replace(r"(L|R)$", "", regex=True)
  557. # For FcMi & FcMa → both-hemisphere bundles → assign pseudo-hemi
  558. df.loc[df["BundleBase"].isin(["FcMi", "FcMa"]), "Hemisphere"] = "B"
  561. # --- DEFINE LAYOUT AND STYLE --------------------------------------------------
  562. tractPosHorz = {
  563. 'AF': (0,0), 'ATR': (0,1), 'CC': (0,2), 'CS': (0,3),
  564. 'IFOF': (1,0), 'ILF': (1,1), 'MLF': (1,2), 'OR': (1,3),
  565. 'SLF': (2,0), 'UNC': (2,1), 'VOF': (2,2), 'pAF': (2,3),
  566. 'FcMi': (3,1), 'FcMa': (3,2)
  567. }
  569. tracts = ['AF', 'ATR', 'CC', 'CS', 'IFOF', 'ILF', 'MLF',
  570. 'OR', 'SLF', 'UNC', 'VOF', 'pAF', 'FcMi', 'FcMa']
  572. colors = ['cyan', 'blue', 'green', 'orange', 'purple', 'brown', 'olive',
  573. 'coral', 'fuchsia', 'yellow', 'indigo', 'violet', 'salmon', 'red']
  575. sns.set(font_scale=2)
  576. sns.set_style("white")
  578. fig, axes = plt.subplots(4, 4, figsize=(14,17), frameon=False)
  579. fig.subplots_adjust(wspace=0.1, hspace=0.9)
  581. import matplotlib.colors as mcolors
  583. def lighten_color(color, amount=0.5):
  584. try:
  585. c = mcolors.cnames[color]
  586. except KeyError:
  587. c = color
  588. rgb = mcolors.to_rgb(c)
  589. return tuple(1 - (1 - x) * amount for x in rgb)
  592. # --- MAIN LOOP ---------------------------------------------------------------
  593. for ct, (bundle, color) in enumerate(zip(tracts, colors)):
  594. if bundle not in tractPosHorz:
  595. continue
  597. ax = axes[tractPosHorz[bundle]]
  599. # Filter data for the base bundle (both hemispheres)
  600. SlopeTract = df.query("BundleBase == @bundle")
  602. # Aggregate per subject & hemisphere
  603. SlopeTractGM = SlopeTract.groupby(
  604. ['subjectID', 'scan_age', 'BundleBase', 'Hemisphere'],
  605. as_index=False, sort=False
  606. )['WeAvGMT1wT2w'].mean()
  608. SlopeTractWM = SlopeTract.groupby(
  609. ['subjectID', 'scan_age', 'BundleBase', 'Hemisphere'],
  610. as_index=False, sort=False
  611. )['t1wt2w'].mean()
  613. # Special case: FcMi & FcMa (ONLY ONE "hemisphere")
  614. if bundle in ["FcMi", "FcMa"]:
  615. sns.scatterplot(
  616. data=SlopeTractWM, x='scan_age', y='t1wt2w',
  617. color=color, s=60, alpha=0.5, ax=ax
  618. )
  619. sns.regplot(
  620. data=SlopeTractWM,
  621. x='scan_age', y='t1wt2w', scatter=False, ax=ax,
  622. line_kws={'color': color, 'alpha': 0.9}
  623. )
  625. sns.scatterplot(
  626. data=SlopeTractGM, x='scan_age', y='WeAvGMT1wT2w',
  627. color='gray', s=60, alpha=0.5, ax=ax
  628. )
  629. sns.regplot(
  630. data=SlopeTractGM,
  631. x='scan_age', y='WeAvGMT1wT2w', scatter=False, ax=ax,
  632. line_kws={'color': 'gray', 'alpha': 0.9}
  633. )
  635. else:
  636. # Normal L/R hemisphere bundles --------------------------------------
  637. light_color = lighten_color(color, amount=0.5)
  639. # WM scatter
  640. sns.scatterplot(
  641. data=SlopeTractWM, x='scan_age', y='t1wt2w',
  642. hue='Hemisphere', style='Hemisphere',
  643. markers={'L': 'o', 'R': 'o'}, s=50,
  644. palette={'L': color, 'R': light_color}, ax=ax,
  645. alpha=0.5, legend=False
  646. )
  648. # WM regression lines
  649. for hemi, line_color in zip(['L', 'R'], [color, light_color]):
  650. sns.regplot(
  651. data=SlopeTractWM[SlopeTractWM['Hemisphere'] == hemi],
  652. x='scan_age', y='t1wt2w', scatter=False, ax=ax,
  653. line_kws={'color': line_color, 'alpha': 0.8}
  654. )
  656. # GM scatter
  657. sns.scatterplot(
  658. data=SlopeTractGM, x='scan_age', y='WeAvGMT1wT2w',
  659. hue='Hemisphere', style='Hemisphere',
  660. markers={'L': 'o', 'R': 'o'}, s=50,
  661. palette={'L': 'darkgray', 'R': 'lightgray'}, ax=ax,
  662. alpha=0.5, legend=False
  663. )
  665. # GM regression
  666. for hemi, line_color in zip(['L', 'R'], ['darkgray', 'lightgray']):
  667. sns.regplot(
  668. data=SlopeTractGM[SlopeTractGM['Hemisphere'] == hemi],
  669. x='scan_age', y='WeAvGMT1wT2w', scatter=False, ax=ax,
  670. line_kws={'color': line_color, 'alpha': 0.8}
  671. )
  674. # --- Titles & labels -----------------------------------------------------
  675. ax.set_title(bundle, pad=5, fontsize=25)
  676. ax.set_xlabel('Scan Age', labelpad=15, fontsize=20)
  678. if ct in [0, 4, 8, 12]: # first column
  679. ax.set_ylabel("T1w/T2w", fontsize=20)
  680. ax.set_xticks([30, 34, 38, 42])
  681. ax.spines['right'].set_visible(False)
  682. ax.spines['top'].set_visible(False)
  684. else:
  685. ax.set_ylabel('T1w/T2w', fontsize=14)
  686. ax.set_xticks([30, 34, 38, 42])
  687. ax.yaxis.set_visible(False)
  688. ax.spines['right'].set_visible(False)
  689. ax.spines['left'].set_visible(False)
  690. ax.spines['top'].set_visible(False)
  692. ax.spines['right'].set_visible(False)
  693. ax.spines['top'].set_visible(False)
  696. # Turn off unused axes
  697. axes[3,3].axis("off")
  698. axes[3,0].axis("off")
  700. # --- SAVE & SHOW -------------------------------------------------------------
  701. fig.savefig('./figures/scanage_per_bundle_bothhemis.png', dpi=600)
  702. plt.show()
  704. # %% [markdown]
  705. # ## R1
  707. # %%
  708. #Compute slope for scan age for R1 of white matter
  710. tractCount=SlopeDataframe['tractID'].unique()
  712. coeff=[] #Creating arrays
  713. se = []
  714. coeffType=[]
  715. tractTrack1=[]
  716. intercept=[]
  717. r2 = []
  718. adjustR2 = []
  719. aic = []
  720. pvals=[]
  722. for tract in tractCount:
  724. dfForStats=SlopeDataframe[(SlopeDataframe['tractID']==tract)]
  725. dfForStats.reset_index(level=0, inplace=True)
  727. scanAgeArray=np.array(dfForStats['age'])
  728. scanAgeArray=scanAgeArray.reshape(1,-1)
  729. scanAgeNorm=preprocessing.normalize(scanAgeArray,axis=1)
  730. scanAgeNorm=scanAgeNorm.reshape(-1)
  731. dfForStats["scanAgeNorm"]=scanAgeNorm
  733. R1Array=np.array(dfForStats['R1'])
  734. R1Array=R1Array.reshape(1,-1)
  735. R1Norm=preprocessing.normalize(R1Array,axis=1)
  736. R1Norm=R1Norm.reshape(-1)
  737. dfForStats["R1Norm"]=R1Norm
  739. df = dfForStats
  740. md = smf.ols(formula="R1~ 1 + age", data=df) #Predicting dti_mdNorm through the other variables
  741. mdf = md.fit() # Actually commanding the linear analysis, it could be any other analysis method
  742. results=mdf.summary()
  744. # Note that tables is a list. The table at index 1 is the "core" table. Additionally, read_html puts dfs in a list, so we want index 0
  745. results_as_html = results.tables[1].as_html()
  746. LM_results=pd.read_html(results_as_html, header=0, index_col=0)[0]
  748. results_as_html_R2 = results.tables[0].as_html()
  749. LM_results_R2=pd.read_html(results_as_html_R2, index_col=0)[0]
  751. results_R2s=LM_results_R2[3]
  752. #LM_results.head()
  753. #LM_results_R2.head()
  754. #results_R2s.head()
  756. SE=LM_results['std err']
  758. inter=mdf.params['Intercept']
  759. intercept.append(inter)
  761. #SE.head()
  762. se.append(SE["age"])
  763. coeff.append(mdf.params['age'])
  764. coeff_type=3
  765. tractTrack1.append(tract)
  766. #nodeTrack1.append(node)
  767. coeffType.append(coeff_type) #Add another line below
  768. r2.append(results_R2s['Dep. Variable:'])
  769. adjustR2.append(results_R2s['Model:'])
  770. aic.append(results_R2s['No. Observations:'])
  771. pvals.append(results_R2s['Date:'])
  773. print(mdf.summary())
  774. df=pd.DataFrame(tractTrack1)
  775. df.columns=['tractID']
  776. df.insert(1,"coeffType",coeffType)
  777. df.insert(2,"coeff",coeff)
  778. df.insert(3,"se",se)
  779. df.insert(4,"intercept",intercept)
  780. df.insert(5,"r2",r2)
  781. df.insert(6,"adjustR2",adjustR2)
  782. df.insert(7,"aic",aic)
  783. df.insert(8,"pvals",pvals)
  785. os.makedirs('./outputData',exist_ok=True)
  786. df.to_csv('./outputData/LM_result_scanAge_noNorm_WM_AverageSubj_R1.csv')
  787. df.head()
  788. df.min(axis=0, numeric_only=True)
  789. df.max(axis=0, numeric_only=True)
  790. df.mean(axis=0, numeric_only=True)
  791. dfScanAge=df
  792. df.head()
  794. statsPerBundle=df.groupby('tractID').mean()
  795. statsPerBundle.head()
  796. statsPerBundle.to_csv('./outputData/LM_result_scanAge_noNorm_meanPerBundle_WM_AverageSubj_R1.csv')
  798. # %%
  799. #Compute slope for scan age for R1 of gray matter
  801. tractCount=SlopeDataframe['tractID'].unique()
  803. coeff=[] #Creating arrays
  804. se = []
  805. coeffType=[]
  806. tractTrack1=[]
  807. intercept=[]
  808. r2 = []
  809. adjustR2 = []
  810. aic = []
  811. pvals=[]
  813. for tract in tractCount:
  815. dfForStats=SlopeDataframe[(SlopeDataframe['tractID']==tract)]
  816. dfForStats.reset_index(level=0, inplace=True)
  818. scanAgeArray=np.array(dfForStats['age'])
  819. scanAgeArray=scanAgeArray.reshape(1,-1)
  820. scanAgeNorm=preprocessing.normalize(scanAgeArray,axis=1)
  821. scanAgeNorm=scanAgeNorm.reshape(-1)
  822. dfForStats["scanAgeNorm"]=scanAgeNorm
  824. R1Array=np.array(dfForStats['GM_R1'])
  825. R1Array=R1Array.reshape(1,-1)
  826. R1Norm=preprocessing.normalize(R1Array,axis=1)
  827. R1Norm=R1Norm.reshape(-1)
  828. dfForStats["R1Norm"]=R1Norm
  830. df = dfForStats
  831. md = smf.ols(formula="GM_R1~ 1 + age", data=df) #Predicting dti_mdNorm through the other variables
  832. mdf = md.fit() # Actually commanding the linear analysis, it could be any other analysis method
  833. results=mdf.summary()
  835. # Note that tables is a list. The table at index 1 is the "core" table. Additionally, read_html puts dfs in a list, so we want index 0
  836. results_as_html = results.tables[1].as_html()
  837. LM_results=pd.read_html(results_as_html, header=0, index_col=0)[0]
  839. results_as_html_R2 = results.tables[0].as_html()
  840. LM_results_R2=pd.read_html(results_as_html_R2, index_col=0)[0]
  842. results_R2s=LM_results_R2[3]
  843. #LM_results.head()
  844. #LM_results_R2.head()
  845. #results_R2s.head()
  847. SE=LM_results['std err']
  849. inter=mdf.params['Intercept']
  850. intercept.append(inter)
  852. #SE.head()
  853. se.append(SE["age"])
  854. coeff.append(mdf.params['age'])
  855. coeff_type=3
  856. tractTrack1.append(tract)
  857. #nodeTrack1.append(node)
  858. coeffType.append(coeff_type) #Add another line below
  859. r2.append(results_R2s['Dep. Variable:'])
  860. adjustR2.append(results_R2s['Model:'])
  861. aic.append(results_R2s['No. Observations:'])
  862. pvals.append(results_R2s['Date:'])
  864. print(mdf.summary())
  865. df=pd.DataFrame(tractTrack1)
  866. df.columns=['tractID']
  867. df.insert(1,"coeffType",coeffType)
  868. df.insert(2,"coeff_GM",coeff)
  869. df.insert(3,"se",se)
  870. df.insert(4,"intercept",intercept)
  871. df.insert(5,"r2",r2)
  872. df.insert(6,"adjustR2",adjustR2)
  873. df.insert(7,"aic",aic)
  874. df.insert(8,"pvals",pvals)
  876. os.makedirs('./outputData',exist_ok=True)
  877. df.to_csv('./outputData/LM_result_scanAge_noNorm_GM_AverageSubj_R1.csv')
  878. df.head()
  879. df.min(axis=0, numeric_only=True)
  880. df.max(axis=0, numeric_only=True)
  881. df.mean(axis=0, numeric_only=True)
  882. dfScanAge=df
  883. df.head()
  885. statsPerBundle=df.groupby('tractID').mean()
  886. statsPerBundle.head()
  887. statsPerBundle.to_csv('./outputData/LM_result_scanAge_noNorm_meanPerBundle_GM_AverageSubj_R1.csv')
  889. # %%
  890. SlopeWM_R1 = pd.read_csv('./outputData/LM_result_scanAge_noNorm_WM_AverageSubj_R1.csv', index_col=0)
  891. SlopeGM_R1 = pd.read_csv('./outputData/LM_result_scanAge_noNorm_GM_AverageSubj_R1.csv', index_col=0)
  893. SlopeBothR1 = pd.merge(SlopeWM_R1, SlopeGM_R1[['tractID', 'coeff_GM', 'r2', 'aic']], on='tractID')
  894. SlopeBothR1
  896. SlopeBothR1=SlopeBothR1.set_index(['tractID']).reindex(['AFL', 'AFR', 'ATRL', 'ATRR', 'CCL', 'CCR', 'CSL', 'CSR', 'FcMi', 'FcMa', 'IFOFL', 'IFOFR', 'ILFL', 'ILFR', 'MLFL', 'MLFR', 'ORL', 'ORR', 'SLFL', 'SLFR', 'UNCL', 'UNCR', 'VOFL', 'VOFR', 'pAFL', 'pAFR']).reset_index()
  897. SlopeBothR1
  899. # %%
  900. # Paired t-test between white matter and gray matter coefficients
  901. t_stat, p_val = stats.ttest_rel(SlopeBothR1['coeff'], SlopeBothR1['coeff_GM'])
  902. print(f"CoeffWM vs. CoeffGM t = {t_stat:.4f}, p = {p_val:.4f}")
  904. # Print mean values
  905. mean_wm = SlopeBothR1['coeff'].mean()
  906. mean_gm = SlopeBothR1['coeff_GM'].mean()
  907. print(f"Mean Coeff WM: {mean_wm:.4f}")
  908. print(f"Mean Coeff GM: {mean_gm:.4f}")
  910. # Calculate difference metrics y=GM, x=WM
  911. SlopeBothR1['diff_r2'] = SlopeBothR1['r2_x'] - SlopeBothR1['r2_y']
  912. SlopeBothR1['diff_aic'] = SlopeBothR1['aic_x'] - SlopeBothR1['aic_y']
  914. # Print mean differences
  915. mean_diff_r2 = SlopeBothR1['diff_r2'].mean()
  916. mean_diff_aic = SlopeBothR1['diff_aic'].mean()
  917. print(f"Mean Δr²: {mean_diff_r2:.4f}")
  918. print(f"Mean ΔAIC: {mean_diff_aic:.4f}")
  920. # %%
  921. sns.set_style('white');
  922. fig1=sns.lmplot(
  923. data=SlopeBothR1, x='coeff_GM', y="coeff", hue="tractID", height=10, scatter_kws={"s": 400}, fit_reg=False, legend=False, palette=color_list_chosen
  924. )
  925. fig1=sns.regplot(data=SlopeBothR1, x='coeff_GM', y="coeff", scatter=False, ax=fig1.axes[0, 0], line_kws={"color": "darkgrey"})
  926. res=scipy.stats.pearsonr(SlopeBothR1['coeff_GM'], SlopeBothR1['coeff'])
  927. print(res)
  928. res.confidence_interval()
  929. fig1.legend(title='Tract', fontsize='15', title_fontsize='20', bbox_to_anchor=(0.7, 0.25, 0.5, 0.5), loc='right', borderaxespad=0, frameon=False)
  930. plt.xlabel("R1 Slope in GM [s$^{-1}$/w]", fontsize=50)
  931. plt.ylabel("R1 Slope in WM [s$^{-1}$/w]", fontsize=50)
  932. plt.xticks([0.001, 0.003, 0.005], fontsize=45)
  933. plt.yticks(fontsize=45)
  934. plt.savefig('./figures/R1Slope.png', dpi=600, bbox_inches='tight', pad_inches=0.1)
  936. # %%
  937. import os
  938. import numpy as np
  939. import pandas as pd
  940. import seaborn as sns
  941. import matplotlib.pyplot as plt
  942. from matplotlib.lines import Line2D
  943. import matplotlib.colors as mcolors
  945. # --- PREPARE DATA -------------------------------------------------------------
  946. df = SlopeDataframe.copy()
  948. # Extract Hemisphere from tract names (e.g. "AF_L" → "L"); midline bundles will become NaN here
  949. df["Hemisphere"] = df["tractID"].str.extract(r'(L|R)$')
  951. # Extract base bundle name (e.g. "AF_L" → "AF")
  952. df["BundleBase"] = df["tractID"].str.replace(r'(L|R)$', '', regex=True)
  954. # --- DEFINE LAYOUT AND STYLE --------------------------------------------------
  955. tractPosHorz = {
  956. 'AF': (0,0), 'ATR': (0,1), 'CC': (0,2), 'CS': (0,3),
  957. 'IFOF':(1,0), 'ILF':(1,1), 'MLF':(1,2), 'OR':(1,3),
  958. 'SLF':(2,0), 'UNC':(2,1), 'VOF':(2,2), 'pAF':(2,3),
  959. 'FcMi':(3,1), 'FcMa':(3,2)
  960. }
  962. tracts = ['AF', 'ATR', 'CC', 'CS', 'IFOF', 'ILF', 'MLF',
  963. 'OR', 'SLF', 'UNC', 'VOF', 'pAF', 'FcMi', 'FcMa']
  965. colors = ['cyan', 'blue', 'green', 'orange', 'purple', 'brown', 'olive',
  966. 'coral', 'fuchsia', 'yellow', 'indigo', 'violet', 'salmon', 'red']
  968. sns.set(font_scale=2)
  969. sns.set_style("white")
  971. fig, axes = plt.subplots(4, 4, figsize=(14,17), frameon=False)
  972. fig.subplots_adjust(wspace=0.1, hspace=0.9)
  974. def lighten_color(color, amount=0.5):
  975. """
  976. Lightens the given color by blending it with white.
  977. amount=0 → white, amount=1 → original color.
  978. """
  979. try:
  980. c = mcolors.cnames[color]
  981. except KeyError:
  982. c = color
  983. rgb = mcolors.to_rgb(c)
  984. return tuple(1 - (1 - x) * amount for x in rgb)
  986. # --- MAIN LOOP ---------------------------------------------------------------
  987. midline_bundles = ["FcMi", "FcMa"] # bundles without L/R hemispheres
  989. for ct, (bundle, color) in enumerate(zip(tracts, colors)):
  990. # skip if not in layout mapping
  991. if bundle not in tractPosHorz:
  992. continue
  994. ax = axes[tractPosHorz[bundle]]
  996. # Filter data for this bundle base (works for both midline and hemispheric bundles)
  997. SlopeTract = df.query("BundleBase == @bundle")
  999. # If there's no data at all for this bundle, skip plotting (avoid empty panels)
  1000. if SlopeTract.shape[0] == 0:
  1001. ax.text(0.5, 0.5, "No data", ha='center', va='center', transform=ax.transAxes, fontsize=14)
  1002. ax.set_title(bundle, pad=5, fontsize=20)
  1003. ax.set_xticks([])
  1004. ax.set_yticks([])
  1005. continue
  1007. # --- Different grouping for midline vs. hemispheric bundles ----------------
  1008. if bundle in midline_bundles:
  1009. # Group WITHOUT 'Hemisphere' because it's NaN for midline bundles
  1010. SlopeTractWM = (
  1011. SlopeTract
  1012. .groupby(['subjectID', 'age', 'BundleBase'], as_index=False, sort=False)['R1']
  1013. .mean()
  1014. )
  1015. SlopeTractGM = (
  1016. SlopeTract
  1017. .groupby(['subjectID', 'age', 'BundleBase'], as_index=False, sort=False)['GM_R1']
  1018. .mean()
  1019. )
  1021. # Mark as midline for plotting convenience
  1022. SlopeTractWM["Hemisphere"] = "M"
  1023. SlopeTractGM["Hemisphere"] = "M"
  1025. # Plot WM (single color)
  1026. sns.scatterplot(
  1027. data=SlopeTractWM, x='age', y='R1',
  1028. color=color, s=50, alpha=0.5, ax=ax
  1029. )
  1030. sns.regplot(
  1031. data=SlopeTractWM, x='age', y='R1',
  1032. scatter=False, ax=ax,
  1033. line_kws={'color': color, 'alpha': 0.8}
  1034. )
  1036. # Plot GM (single gray color)
  1037. sns.scatterplot(
  1038. data=SlopeTractGM, x='age', y='GM_R1',
  1039. color='gray', s=50, alpha=0.5, ax=ax
  1040. )
  1041. sns.regplot(
  1042. data=SlopeTractGM, x='age', y='GM_R1',
  1043. scatter=False, ax=ax,
  1044. line_kws={'color': 'gray', 'alpha': 0.8}
  1045. )
  1047. else:
  1048. # Hemispheric bundles: group including Hemisphere field
  1049. # Note: groupby will ignore NaN Hemispheres—this is fine because hemispheric bundles should have L/R
  1050. SlopeTractWM = SlopeTract.groupby(
  1051. ['subjectID', 'age', 'BundleBase', 'Hemisphere'],
  1052. as_index=False, sort=False
  1053. )['R1'].mean()
  1055. SlopeTractGM = SlopeTract.groupby(
  1056. ['subjectID', 'age', 'BundleBase', 'Hemisphere'],
  1057. as_index=False, sort=False
  1058. )['GM_R1'].mean()
  1060. light_color = lighten_color(color, amount=0.5)
  1062. # White matter scatter (L/R)
  1063. sns.scatterplot(
  1064. data=SlopeTractWM, x='age', y='R1', hue='Hemisphere', style='Hemisphere',
  1065. markers={'L': 'o', 'R': 'o'}, s=50,
  1066. palette={'L': color, 'R': light_color}, ax=ax,
  1067. alpha=0.5, legend=False
  1068. )
  1070. # Regression lines for L and R (if any data present)
  1071. for hemi, line_color in zip(['L','R'], [color, light_color]):
  1072. hemi_df = SlopeTractWM[SlopeTractWM['Hemisphere'] == hemi]
  1073. if hemi_df.shape[0] > 1: # need at least 2 points for a regression line
  1074. sns.regplot(
  1075. data=hemi_df, x='age', y='R1', scatter=False, ax=ax,
  1076. line_kws={'color': line_color, 'alpha': 0.8}
  1077. )
  1079. # Gray matter scatter & lines
  1080. sns.scatterplot(
  1081. data=SlopeTractGM, x='age', y='GM_R1', hue='Hemisphere', style='Hemisphere',
  1082. markers={'L': 'o', 'R': 'o'}, s=50,
  1083. palette={'L':'darkgray', 'R':'lightgray'}, ax=ax,
  1084. alpha=0.5, legend=False
  1085. )
  1087. for hemi, line_color in zip(['L','R'], ['darkgray','lightgray']):
  1088. hemi_gm = SlopeTractGM[SlopeTractGM['Hemisphere'] == hemi]
  1089. if hemi_gm.shape[0] > 1:
  1090. sns.regplot(
  1091. data=hemi_gm, x='age', y='GM_R1', scatter=False, ax=ax,
  1092. line_kws={'color': line_color, 'alpha': 0.8}
  1093. )
  1095. # --- TITLES AND LABELS -----------------------------------------------------
  1096. ax.set_title(bundle, pad=5, fontsize=25)
  1097. ax.set_xlabel('Scan Age', labelpad=15, fontsize=20)
  1099. if ct in [0, 4, 8, 12]: # first column
  1100. ax.set_ylabel("R1 [s$^{-1}$]", fontsize=20)
  1101. ax.set_xticks([40, 42, 44, 46])
  1102. ax.set_xticklabels([40, 42, 44, 46], fontsize=20)
  1103. ax.set_yticks([0.4, 0.45, 0.5])
  1104. ax.set_yticklabels([0.4, 0.45, 0.5], fontsize=20)
  1105. ax.spines['right'].set_visible(False)
  1106. ax.spines['top'].set_visible(False)
  1107. else:
  1108. ax.set_ylabel('R1 [s$^{-1}$]', fontsize=14)
  1109. ax.set_xticks([40, 42, 44, 46])
  1110. ax.set_xticklabels([40, 42, 44, 46], fontsize=20)
  1111. ax.set_yticks([0.4, 0.45, 0.5])
  1112. ax.yaxis.set_visible(False)
  1113. ax.spines['right'].set_visible(False)
  1114. ax.spines['left'].set_visible(False)
  1115. ax.spines['top'].set_visible(False)
  1117. # Clean frame (ensure spines hidden where desired)
  1118. ax.spines['right'].set_visible(False)
  1119. ax.spines['top'].set_visible(False)
  1121. # Turn off unused axes
  1122. axes[3,3].axis("off")
  1123. axes[3,0].axis("off")
  1125. # --- SAVE & SHOW -------------------------------------------------------------
  1126. os.makedirs('./figures', exist_ok=True)
  1127. plt.savefig('./figures/R1ScanAge_Hemispheres.png', dpi=600, bbox_inches='tight')
  1128. plt.show()

Figure3.ipynb at commit 57e9055, no license · at the source

Overview

Authors: Stephanie Zika1,2,3, Kelly Chang4,5, Altan Orhon4,5, John Kruper4,5, Christina Tyagi6, Xiaoqian Yan7, Sarah Tung6, Kalanit Grill-Spector6,8, Ariel Rokem4,5, Mareike Grotheer1,2,3
ORCID iDs: Stephanie Zika, Kelly Chang, John Kruper, Christina Tyagi, Xiaoqian Yan, Sarah Tung, Kalanit Grill-Spector, Ariel Rokem, Mareike Grotheer
  1. Department of Psychology, Philipps-Universität Marburg, Marburg, Germany
  2. Center for Mind, Brain and Behavior - CMBB, Universities of Marburg, Gießen, Darmstadt, Germany
  3. Lise Meitner Research Group Neuroplasticity in Development and Learning, Max Planck Institute for Human Development, Berlin, Germany
  4. Department of Psychology, University of Washington, Seattle, WA USA
  5. eScience Institute, University of Washington, Seattle, WA USA
  6. Department of Psychology, Stanford University, Stanford, CA USA
  7. Institute of Science and Technology for Brain-Inspired Intelligence, Fudan University, Shanghai, China
  8. Wu Tsai Neurosciences Institute, Stanford University, Stanford, CA USA
Institutions: Philipps University of Marburg (Germany); Max Planck Institute for Human Development (Germany); University of Washington (United States); Stanford University (United States); Shanghai Medical College of Fudan University (China); Fudan University (China)
Journal: Nature communications, volume 17, issue 1, article 5353
Dates: received 8 July 2025; accepted 7 May 2026; published online 17 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-73366-9 · PMID 42310317 · PMCID PMC13276173 · OpenAlex W7165007514
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), developmental (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Preprocessing, Connectivity, fMRI & imaging
Keywords: Human behaviour, Developmental neurogenesis, Neuronal development
MeSH: Cerebral Cortex*, Gray Matter*, Myelin Sheath*, White Matter*, Child, Preschool, Female, Humans, Infant, Infant, Newborn, Magnetic Resonance Imaging, Male (* major topic)
Topic: Advanced Neuroimaging Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: NIH National Institutes of Health & National Science Foundation (R01-EY033835, MH121868, R01-EB027585, R01-EY033628); European Research Council (319456); NIH NatioNIH National Institutes of Health & National Science Foundationnal Institutes of Health & National Science Foundation (MH121867)
Citations: cited by 1 paper (Europe PMC); 131 references in the paper

Abstract

The infant brain undergoes rapid myelination that is critical for healthy brain function. This development has been characterized for gray and white matter independently, but the link between gray and white matter myelination remains unexplored. To close this knowledge gap, we evaluated two complementary myelin-sensitive imaging metrics: Large-scale (N = 273) T1w/T2w and quantitative (N = 21) R1 data. Automated software was employed to identify 26 white matter bundles and map their cortical terminations, before evaluating T1w/T2w and R1 development shortly after birth. Here we show that for both metrics mean values as well as developmental slopes are correlated across tissues. The synchrony of brain T1w/T2w is impacted by postmenstrual age and prematurity, whereas inter-individual differences in this synchrony predict motor outcomes at 17 − 25 months of age. As T1w/T2w and R1 are associated with myelin content, our results reveal an intricate relationship between gray and white matter myelination.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

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

EduNeuroLab/WMGMMyelinInfants

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 57e9055dc133eca1fdbe06690a8fd9d233fb5e5c, 31 March 2026
Languages: Jupyter (8), Python (1)
Size: 12 files, 9 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 8 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: pandas (6 files), Matplotlib (5 files), NumPy (5 files), SciPy (5 files), seaborn (5 files), NiBabel (4 files), FreeSurfer (3 files), statsmodels (3 files), MRtrix3 (2 files), Nipype (2 files), scikit-learn (2 files), DIPY (1 file), scikit-image (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
10 files
At the source: github.com/EduNeuroLab/WMGMMyelinInfants

Code availability

All data were analyzed using open-source software, including MRtrix110,128 and pyBabyAFQ, which we shared as a component of pyAFQ (https://tractometry.org/pyAFQ)124,129,130. Code that implements the tractography pipeline, example code to perform bundle identification with pyBabyAFQ, and code used to generate the main figures of this manuscript are made available in GitHub (https://github.com/EduNeuroLab/WMGMMyelinInfants)131.

Reproduced under the paper's license (CC BY), from the paper cited above.

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;
  • 9 scripts, each with its path and the digest of its content;
  • 11 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

Datasets cited

Data Availability Statement

All data required to generate the main figures are provided as a Source Data file with this paper and are also made available in GitHub (https://github.com/EduNeuroLab/WMGMMyelinInfants). In addition, data for analyses within this work are provided in Figshare (https://figshare.com/articles/dataset/Data_from_Cortical_and_white_matter_myelination_proceed_in_concert_during_early_infancy_/30667400), in accordance with dHCP data use terms and conditions. Source data are provided with this paper.

All data were analyzed using open-source software, including MRtrix110,128 and pyBabyAFQ, which we shared as a component of pyAFQ (https://tractometry.org/pyAFQ)124,129,130. Code that implements the tractography pipeline, example code to perform bundle identification with pyBabyAFQ, and code used to generate the main figures of this manuscript are made available in GitHub (https://github.com/EduNeuroLab/WMGMMyelinInfants)131.

Reproduced under the paper's license (CC BY), from the paper cited above.

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, 10 authors, 3 keywords, 11 MeSH terms, 3 funders, 122 references.

Cite

This paper

Zika, S., Chang, K., Orhon, A., Kruper, J., Tyagi, C., Yan, X., Tung, S., Grill-Spector, K., Rokem, A., & Grotheer, M. (2026). Cortical and white matter myelination proceed in concert during early infancy. Nature communications, 17(1), 5353. https://doi.org/10.1038/s41467-026-73366-9

BibTeX

@article{zika2026cortical,
author = {Zika, Stephanie and Chang, Kelly and Orhon, Altan and Kruper, John and Tyagi, Christina and Yan, Xiaoqian and Tung, Sarah and Grill-Spector, Kalanit and Rokem, Ariel and Grotheer, Mareike},
title = {{Cortical and white matter myelination proceed in concert during early infancy}},
journal = {Nature communications},
year = {2026},
month = jun,
volume = {17},
number = {1},
pages = {5353},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-73366-9},
url = {https://doi.org/10.1038/s41467-026-73366-9},
pmid = {42310317},
pmcid = {PMC13276173}
}

RIS

TY - JOUR
AU - Zika, Stephanie
AU - Chang, Kelly
AU - Orhon, Altan
AU - Kruper, John
AU - Tyagi, Christina
AU - Yan, Xiaoqian
AU - Tung, Sarah
AU - Grill-Spector, Kalanit
AU - Rokem, Ariel
AU - Grotheer, Mareike
TI - Cortical and white matter myelination proceed in concert during early infancy
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/06/17
VL - 17
IS - 1
SP - 5353
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-73366-9
UR - https://doi.org/10.1038/s41467-026-73366-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-73366-9",
"type": "article-journal",
"title": "Cortical and white matter myelination proceed in concert during early infancy",
"container-title": "Nature communications",
"author": [
{
"family": "Zika",
"given": "Stephanie"
},
{
"family": "Chang",
"given": "Kelly"
},
{
"family": "Orhon",
"given": "Altan"
},
{
"family": "Kruper",
"given": "John"
},
{
"family": "Tyagi",
"given": "Christina"
},
{
"family": "Yan",
"given": "Xiaoqian"
},
{
"family": "Tung",
"given": "Sarah"
},
{
"family": "Grill-Spector",
"given": "Kalanit"
},
{
"family": "Rokem",
"given": "Ariel"
},
{
"family": "Grotheer",
"given": "Mareike"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "5353",
"DOI": "10.1038/s41467-026-73366-9",
"PMID": "42310317",
"PMCID": "PMC13276173",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-73366-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
17
]
]
}
}

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.1126/sciadv.aec2348 [code]
Congenital blindness reduces myelination in human visual cortex.
Journal: Science advances
In common: FreeSurfer, NiBabel, pandas, 1 other tool, 19 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, FreeSurfer, statsmodels, 5 other tools, 11 references
[3] doi:10.1162/imag.a.1341 [code]
Massively parallelized brain tractography using compute clusters, supercomputers, and graphics processing units.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: DIPY, NiBabel, pandas, 3 other tools, 5 references, 2 authors
[4] doi:10.1093/cercor/bhag132 [code]
Spatiotemporal white-matter development across early childhood.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: DIPY, MRtrix3, scikit-image, 5 other tools, developmental, 9 references
[5] doi:10.1038/s41467-026-73072-6 [code]
Mapping the spatiotemporal continuum of structural connectivity development across the human connectome in youth.
Journal: Nature communications
In common: MRtrix3, FreeSurfer, NiBabel, 3 other tools, developmental, 12 references
[6] doi:10.1038/s41597-026-06869-1 [code]
Individual Brain Charting: fifth release of high-resolution fMRI data for cognitive mapping.
Journal: Scientific data
In common: DIPY, MRtrix3, Nipype, 9 other tools, author Kalanit Grill-Spector
[7] doi:10.1016/j.dcn.2026.101775 [code]
Neonatal brain-age models in full- and preterm infants.
Journal: Developmental cognitive neuroscience
In common: NiBabel, statsmodels, seaborn, 5 other tools, developmental, 9 references
[8] doi:10.1038/s41598-026-51531-w [code]
Multimodal age-dependent diffusion-MRI analysis of the neocortex in a rat model of cortical dysplasia.
Journal: Scientific reports
In common: DIPY, MRtrix3, scikit-image, 6 other tools, 6 references
[9] doi:10.1002/hbm.70562 [code]
Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning.
Journal: Human brain mapping
In common: MRtrix3, NiBabel, statsmodels, 6 other tools, 6 references
[10] doi:10.1038/s42003-026-10276-y [code]
The cellular correlates and adolescent reorganisation of cortical myelination networks in the common marmoset.
Journal: Communications biology
In common: Nipype, FreeSurfer, NiBabel, 7 other tools, 5 references

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.