OSCR

Dissecting heterogeneous brain development and aging using voxelwise normative models

Code ↔ Paper

12 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 12 matches
  1. [1] § Methods › Datasets, preprocessing and brain measures extraction ↔ voxelwise_models_estimation/02_prepare_voxelwise_normdata.py, lines 1 to 60 · score 0.79 · HCP Aging, HCP Dev, CMI HBN, IXI, PING, PNC · Report a mismatch
  2. [2] § Methods › Voxelwise normative models › Application to spinocerebellar ataxia type 1 and 3 ↔ voxelwise_models_application/06_KernelRidgeRegression_SCA_patients.py, lines 1 to 40 · score 0.76 · kernel ridge regression, brain voxelwise, ataxia, raw, prediction, spinocerebellar · Report a mismatch
  3. [3] § Methods › Datasets, preprocessing and brain measures extraction ↔ voxelwise_models_estimation/01_preprocess_T1_ants_fsl.py, lines 47 to 101 · score 0.76 · nonlinear transformations, combined affine, Jacobian determinants, geometric, preprocessing, log · Report a mismatch
  4. [4] § Methods › Voxelwise normative models › Application to spinocerebellar ataxia type 1 and 3 ↔ voxelwise_models_application/06_KernelRidgeRegression_SCA_patients.py, lines 1 to 40 · score 0.67 · SCA patients, spinocerebellar ataxia, Mexico, SCA1, SCA3, Nijmegen · Report a mismatch
  5. [5] § Methods › Voxelwise normative models › Model estimation ↔ voxelwise_models_application/05_preterm_subjects_comparison.py, lines 1 to 40 · score 0.66 · brain burden, controls born, Mann, TFCE, extreme, normative models · Report a mismatch
  6. [6] § Methods › Datasets, preprocessing and brain measures extraction ↔ voxelwise_models_estimation/02_prepare_voxelwise_normdata.py, lines 1 to 60 · score 0.64 · CMI HBN, MNI152NLin2009cSym, brain extraction, preprocessed, ANTs, ABCD · Report a mismatch
  7. [7] § Methods › Datasets, preprocessing and brain measures extraction ↔ voxelwise_models_estimation/01_preprocess_T1_ants_fsl.py, lines 123 to 148 · score 0.64 · segmented grey, FSL, FAST, transforms, anat, ANTs · Report a mismatch
  8. [8] § Methods › Voxelwise normative models › Application to spinocerebellar ataxia type 1 and 3 ↔ voxelwise_models_application/06_KernelRidgeRegression_SCA_patients.py, lines 311 to 397 · score 0.63 · age residualized, predictive pipeline, Mexico, Nijmegen, regressor, raw · Report a mismatch
  9. [9] § Results › Preterm birth gives rise to widespread and persistent patterns of deviation ↔ voxelwise_models_application/05_preterm_subjects_comparison.py, lines 1 to 40 · score 0.60 · brain burden, controls born, born preterm, extreme, normative model, scores · Report a mismatch
  10. [10] § Methods › Datasets, preprocessing and brain measures extraction ↔ voxelwise_models_estimation/02_prepare_voxelwise_normdata.py, lines 189 to 272 · score 0.57 · FSL, brain extraction, failures, segmentation, anat, ANTs · Report a mismatch
  11. [11] § Results › Preterm birth gives rise to widespread and persistent patterns of deviation ↔ voxelwise_models_estimation/04_evaluate_voxelwise_models.py, lines 212 to 219 · score 0.53 · internal capsule, cerebellar white matter, voxels · Report a mismatch
  12. [12] § Results › SCA individual-level deviation maps show distinct heterogeneity ↔ voxelwise_models_application/06_KernelRidgeRegression_SCA_patients.py, lines 199 to 259 · score 0.51 · kernel ridge regression, Mexico, SCA1, SCA3, Nijmegen, raw · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 542 lines · 21 KB · no license · 4 matches

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. Created on Mon Aug 19 16:35:32 2024
  6. @author: Alice Chavanne
  8. This script conducts a kernel ridge regression analysis evaluated out-of-sample in patients with spinocerebellar ataxia type 1 and 3,
  9. from two independent datasets. The predictive performance of raw, whole-brain voxelwise measures (log jacobian, GMV & WMV combined ) and normative measures
  10. (Z-scores) are explored separately using the options set.
  12. """
  14. import numpy as np
  15. from pcntoolkit.dataio.fileio import save as ptksave
  16. import os
  17. import pandas as pd
  18. import pickle
  19. from sklearn.preprocessing import StandardScaler, FunctionTransformer
  20. from sklearn.model_selection import LeaveOneOut, GridSearchCV
  21. from sklearn.metrics import r2_score, root_mean_squared_error
  22. from sklearn.pipeline import Pipeline
  23. from sklearn.kernel_ridge import KernelRidge
  24. from sklearn.base import BaseEstimator, TransformerMixin
  25. import matplotlib.pyplot as plt
  26. from scipy.stats import pearsonr
  27. import nibabel as nib
  28. from nibabel import processing
  29. import nilearn as nil
  30. from nilearn import image
  31. from pcntoolkit import dataio
  32. import joblib
  33. from joblib import Parallel, delayed
  36. data_dir = '/path/to/data/dir/with/raw/ants/outputs/'
  37. out_dir = os.path.join(data_dir,'clinical_predictions')
  38. extension = ''#optional extra suffix added to output folder/filenames
  39. model_dir = '/path/to/dir/where/normative/models/live'
  40. datasets = ['SCA1_Nijmegen', 'SCA3_Mexico'] #some indices are hardcoded based on this order so have be careful
  43. #%% Choose options for regression
  45. target = {
  46. "BL_SARA" : 0,
  47. "Y1_SARA" : 0,
  48. "Slope_BL_Y1_SARA" : 1,
  49. }
  51. train_and_test = {
  52. "test_Nijmegen_CV":1,
  53. "test_Nijmegen_and_Mexico_CV" : 0,
  54. }
  56. featureset = {
  57. "BL_wholebrain_log_jacs_Z" : 1,#can increase RAM if 1 here
  58. "BL_wholebrain_log_jacs_raw" : 0,
  59. "BL_wholebrain_mod_gmv_and_wmv_Z": 0,
  60. "BL_wholebrain_mod_gmv_and_wmv_raw": 0,
  61. }
  63. options = {
  64. "gridsearch" : 0,
  65. }
  67. #account for features with non-standardized values, those that need smoothing or resample
  68. if featureset["BL_wholebrain_mod_gmv_and_wmv_raw"] :
  69. options["Xstandard"] = 1
  70. resample = True
  71. smooth = True
  72. elif featureset["BL_wholebrain_log_jacs_raw"] :
  73. options["Xstandard"] = 1
  74. resample = True
  75. smooth = False
  76. else :
  77. options["Xstandard"] = 0
  78. resample = False
  79. smooth = False
  81. seed = 42
  83. filename = "_".join([key for key, value in target.items() if value] +
  84. ["KRR_with"]+
  85. [key for key, value in featureset.items() if value] +
  86. [key for key, value in train_and_test.items() if value] +
  87. [key for key, value in options.items() if value]+
  88. [extension])
  90. final_dir = os.path.join(out_dir,filename)
  91. os.makedirs(final_dir, exist_ok=True)
  92. os.chdir(final_dir)
  95. #%% Nifti specifics
  97. if (featureset["BL_wholebrain_log_jacs_raw"] or featureset["BL_wholebrain_log_jacs_Z"]):
  98. mask_nii = os.path.join(model_dir,'tpl-MNI152NLin2009cSym/tpl-MNI152NLin2009cSym_res-2_T1w_mask_BrainExtractionBrain.nii.gz')
  99. ex_nii = mask_nii
  100. proc_dir = os.path.join(model_dir,'models','BLR_aff_nonlin_log_jacs_wholebrain_2mm_14')
  102. if (featureset["BL_wholebrain_mod_gmv_and_wmv_raw"] or featureset["BL_wholebrain_mod_gmv_and_wmv_Z"]): # combine the two masks
  103. mask_nii_gm = os.path.join(model_dir,'tpl-MNI152NLin2009cSym/tpl-MNI152NLin2009cSym_res-2_label-GM_mask_probseg0.5.nii.gz')
  104. mask_nii_wm = os.path.join(model_dir,'tpl-MNI152NLin2009cSym/tpl-MNI152NLin2009cSym_res-2_label-WM_mask_probseg0.5.nii.gz')
  105. mask_nii_combined = os.path.join(model_dir,'tpl-MNI152NLin2009cSym/tpl-MNI152NLin2009cSym_res-2_combined_label-GM_WM_mask.nii.gz')
  106. proc_dir_gm= os.path.join(model_dir,'models','BLR_aff_nonlin_mod_gmv_labelmask_2mm_smoothed_fwhm6_07')
  107. proc_dir_wm= os.path.join(model_dir,'models','BLR_aff_nonlin_mod_wmv_labelmask_2mm_smoothed_fwhm6_07')
  108. proc_dir = proc_dir_gm #just to grab the demographics - gm and wm are identical on that front
  110. #just for plotting back onto the combined mask
  111. ex_nii = mask_nii_combined
  112. mask_nii = mask_nii_combined
  114. #details for brain data extraction from niftis
  115. vox_size = [2, 2, 2] #2mm
  116. fwhm = 6 #if smooth is true
  117. vox_batch_size = 150
  119. #%% Grab X and y data
  121. #utils
  122. def load_file(subpath, mask_nii, resample, smooth, vox_size):
  124. img = nib.load(subpath)
  125. if resample :
  126. img = nib.processing.resample_to_output(img, vox_size)
  127. if smooth:
  128. img = nil.image.smooth_img(img,fwhm)
  129. data = img.get_fdata()
  130. data = dataio.fileio.vol2vec(data, nib.load(mask_nii).get_fdata())
  132. return data
  134. # prepare for gmv and wmv, load GM and WM images, filter them by their respective masks, and combine them into a single vector consistent with the combined mask.
  135. def load_files_and_combine(sub_path_gm, sub_path_wm, mask_nii_gm, mask_nii_wm, mask_nii_combined,resample, smooth, vox_size) :
  137. gm_img = nib.load(sub_path_gm)
  138. wm_img = nib.load(sub_path_wm)
  140. if resample:
  141. gm_img = nib.processing.resample_to_output(gm_img, vox_size)
  142. wm_img = nib.processing.resample_to_output(wm_img, vox_size)
  143. if smooth:
  144. gm_img = image.smooth_img(gm_img, fwhm)
  145. wm_img = image.smooth_img(wm_img, fwhm)
  147. gm_data = gm_img.get_fdata()
  148. wm_data = wm_img.get_fdata()
  150. gm_mask = nib.load(mask_nii_gm).get_fdata().astype(bool)
  151. wm_mask = nib.load(mask_nii_wm).get_fdata().astype(bool)
  152. combined_mask = nib.load(mask_nii_combined).get_fdata().astype(bool)
  154. combined_data = np.zeros_like(gm_data).astype(float)
  155. combined_data[gm_mask] = gm_data[gm_mask]
  156. combined_data[wm_mask] = wm_data[wm_mask]
  158. data = combined_data[combined_mask]
  159. return data
  162. # now grab data from all datasets
  163. df_sca = []
  164. y_sca = []
  165. X_sca = []
  167. for dataset in datasets:
  168. #grab subjects and clinical scores
  169. df_clin = pd.read_csv(os.path.join(proc_dir, 'transfer_' + dataset, 'dem_patients.csv'))
  171. #grab y
  172. if target["BL_SARA"]:
  173. y_clin = df_clin['SARA_score']
  174. if target["Y1_SARA"]:
  175. if 'Y1_SARA' in df_clin.columns:
  176. df_clin.dropna(subset=['Y1_SARA'], inplace=True)
  177. df_clin= df_clin.reset_index()
  178. y_clin = df_clin.get('Y1_SARA', pd.Series([], name='Y1_SARA'))
  179. else :
  180. df_clin = pd.DataFrame()
  181. if target["Slope_BL_Y1_SARA"]:
  182. if 'Slope_BL_Y1_SARA' in df_clin.columns:
  183. df_clin.dropna(subset=['Slope_BL_Y1_SARA'], inplace=True)
  184. df_clin= df_clin.reset_index()
  185. y_clin = df_clin.get('Slope_BL_Y1_SARA', pd.Series([], name='Slope_BL_Y1_SARA'))
  186. else :
  187. df_clin = pd.DataFrame()
  189. sub_list = df_clin['sub_id'] if 'sub_id' in df_clin.columns else pd.Series([], dtype=str)
  191. # ensure this doesn't break even for datasets with missing timepoints
  192. if len(sub_list) == 0:
  193. df_sca.append(df_clin)
  194. y_sca.append(pd.Series([], name='y_empty'))
  195. X_sca.append(pd.DataFrame())
  196. continue
  199. # grab X
  200. if (featureset["BL_wholebrain_log_jacs_raw"]) :
  201. file_name = 'T1_BrainNorm_Jacobian.nii.gz'
  202. if dataset == 'SCA1_Nijmegen' :
  203. f_pattern = os.path.join(data_dir, 'subjects_data','{}','ses-01','anat','T1_MNI152NLin2009cSym', file_name)
  204. if dataset == 'SCA3_Mexico' :
  205. f_pattern =os.path.join(data_dir,'additional_SCA3_data', 'T1','{}','T1_MNI152NLin2009cSym', file_name)
  207. filepaths_clin = [f_pattern.format(sub) for sub in sorted(sub_list.tolist())]
  208. resp_data_clin = Parallel(n_jobs=-1,verbose =5)(delayed(load_file)(sub_path, mask_nii,resample, smooth, vox_size)for sub_path in filepaths_clin)
  209. resp_data_clin = np.stack(resp_data_clin, axis=0)
  211. if (featureset["BL_wholebrain_mod_gmv_and_wmv_raw"]) :
  212. file_name_gm = 'T1_BrainNorm_GMV.nii.gz'
  213. file_name_wm = 'T1_BrainNorm_WMV.nii.gz'
  215. if dataset == 'SCA1_Nijmegen' :
  216. gm_pattern = os.path.join(data_dir, 'subjects_data','{}','ses-01','anat','T1_MNI152NLin2009cSym', file_name_gm)
  217. wm_pattern = os.path.join(data_dir, 'subjects_data','{}','ses-01','anat','T1_MNI152NLin2009cSym', file_name_wm)
  218. if dataset == 'SCA3_Mexico' :
  219. gm_pattern =os.path.join(data_dir,'additional_SCA3_data', 'T1','{}','T1_MNI152NLin2009cSym', file_name_gm)
  220. wm_pattern =os.path.join(data_dir,'additional_SCA3_data', 'T1','{}','T1_MNI152NLin2009cSym', file_name_wm)
  222. resp_data_clin = Parallel(n_jobs=-1,verbose =5)(delayed(load_files_and_combine)(
  223. gm_pattern.format(sub), wm_pattern.format(sub),
  224. mask_nii_gm, mask_nii_wm, mask_nii_combined,
  225. resample, smooth, vox_size)
  226. for sub in sorted(sub_list.tolist()))
  227. resp_data_clin = np.stack(resp_data_clin, axis=0)
  229. if (featureset["BL_wholebrain_log_jacs_Z"]) :
  230. mask = dataio.fileio.load(mask_nii, mask=mask_nii, vol=False).T
  231. n_voxels = len(mask)
  232. n_batches = n_voxels // vox_batch_size + int(n_voxels % vox_batch_size != 0)
  233. voxel_indices = list(range(n_voxels))
  235. Z_all = []
  236. for b in range(n_batches)[:]:
  237. Z_batch = []
  238. batch = f'batch_{b}'
  239. #grab Zs from normdata
  240. with open(os.path.join(proc_dir,'transfer_' + dataset, batch, 'norm_data_clin.pkl'), 'rb') as f:
  241. norm_data= pickle.load(f)
  242. norm_data.load_zscores(save_dir = os.path.join(proc_dir,'transfer_' + dataset, batch,'results'))
  244. for sub in sub_list.tolist() :
  245. idx = np.where(norm_data.subject_ids.values == sub)[0]
  246. sub_batch = norm_data.isel(observations=idx)
  248. Z_sub = pd.DataFrame(sub_batch.Z.values, columns=sub_batch.Z.response_vars.values)
  249. Z_batch.append(Z_sub)
  250. Z_batch = pd.concat(Z_batch)
  251. Z_all.append(Z_batch)
  253. if len(Z_all)>0 :
  254. resp_data_clin = pd.concat(Z_all, axis = 1)
  255. else :
  256. resp_data_clin = Z_all
  259. if (featureset["BL_wholebrain_mod_gmv_and_wmv_Z"]) :
  261. #grab GM Zs from normdata
  262. mask = dataio.fileio.load(mask_nii_gm, mask=mask_nii_gm, vol=False).T
  263. n_voxels = len(mask)
  264. n_batches = n_voxels // vox_batch_size + int(n_voxels % vox_batch_size != 0)
  266. Z_all_gm = []
  267. for b in range(n_batches)[:]:
  268. Z_batch = []
  269. batch = f'batch_{b}'
  271. with open(os.path.join(proc_dir_gm,'transfer_' + dataset, batch, 'norm_data_clin.pkl'), 'rb') as f:
  272. norm_data= pickle.load(f)
  273. norm_data.load_zscores(save_dir = os.path.join(proc_dir_gm,'transfer_' + dataset, batch,'results'))
  275. for sub in sub_list.tolist() :
  276. idx = np.where(norm_data.subject_ids.values == sub)[0]
  277. sub_batch = norm_data.isel(observations=idx)
  279. Z_sub = pd.DataFrame(sub_batch.Z.values, columns=sub_batch.Z.response_vars.values)
  280. Z_batch.append(Z_sub)
  281. Z_batch = pd.concat(Z_batch)
  282. Z_all_gm.append(Z_batch)
  284. Z_all_gm = np.array(pd.concat(Z_all_gm, axis=1))
  286. #grab GM Zs from normdata
  287. mask = dataio.fileio.load(mask_nii_wm, mask=mask_nii_wm, vol=False).T
  288. n_voxels = len(mask)
  289. n_batches = n_voxels // vox_batch_size + int(n_voxels % vox_batch_size != 0)
  291. Z_all_wm = []
  292. for b in range(n_batches)[:]:
  293. Z_batch = []
  294. batch = f'batch_{b}'
  296. with open(os.path.join(proc_dir_wm,'transfer_' + dataset, batch, 'norm_data_clin.pkl'), 'rb') as f:
  297. norm_data= pickle.load(f)
  298. norm_data.load_zscores(save_dir = os.path.join(proc_dir_wm,'transfer_' + dataset, batch,'results'))
  300. for sub in sub_list.tolist() :
  301. idx = np.where(norm_data.subject_ids.values == sub)[0]
  302. sub_batch = norm_data.isel(observations=idx)
  304. Z_sub = pd.DataFrame(sub_batch.Z.values, columns=sub_batch.Z.response_vars.values)
  305. Z_batch.append(Z_sub)
  306. Z_batch = pd.concat(Z_batch)
  307. Z_all_wm.append(Z_batch)
  309. Z_all_wm = np.array(pd.concat(Z_all_wm, axis = 1))
  311. #combine GM and WM Zs with mask handling
  312. mask_combined = nib.load(mask_nii_combined)
  313. mask_combined = mask_combined.get_fdata().astype(bool)
  315. mask_gm = nib.load(mask_nii_gm)
  316. mask_gm = mask_gm.get_fdata().astype(bool)
  318. mask_wm = nib.load(mask_nii_wm)
  319. mask_wm = mask_wm.get_fdata().astype(bool)
  321. Z_all_combined = []
  323. for s in range(len(sub_list.tolist())):
  325. #re-project both Z vectors into nifti space
  326. Z_gm = Z_all_gm[s,:]
  327. Z_wm = Z_all_wm[s,:]
  329. #combine them
  330. Z_combined_vol = np.zeros_like(mask_gm).astype(float)
  331. Z_combined_vol[mask_gm] = Z_gm
  332. Z_combined_vol[mask_wm] = Z_wm
  334. #re-vectorize using combined mask
  335. Z_combined = Z_combined_vol[mask_combined]
  336. Z_all_combined.append(Z_combined)
  338. resp_data_clin = Z_all_combined
  340. #if we age-residualize, append the age column to what will be X, the transformer will take care of it
  341. if (featureset["BL_wholebrain_log_jacs_raw"] or featureset["BL_wholebrain_mod_gmv_and_wmv_raw"]):
  342. resp_data_clin = np.concatenate((resp_data_clin, df_clin["age"].to_numpy().reshape(-1, 1)), 1)
  344. df_sca.append(df_clin)
  345. y_sca.append(y_clin)
  346. X_sca.append(pd.DataFrame(resp_data_clin))
  349. #Now, only get those necessary for the prediction
  351. if train_and_test["test_Nijmegen_CV"]:
  352. X_all = X_sca[0]
  353. X_all = pd.DataFrame(X_all)
  354. y_all = y_sca[0]
  355. df_test = df_sca[0] #we take the df for the end residuals graph
  357. if train_and_test["test_Nijmegen_and_Mexico_CV"]:
  358. X_all = pd.concat([X_sca[0],X_sca[1]])
  359. X_all = pd.DataFrame(X_all)
  360. y_all = pd.concat([y_sca[0], y_sca[1]])
  361. y_all = y_all.reset_index(drop=True) #otehrwsie will throw error
  362. df_test = pd.concat([df_sca[0],df_sca[1]]) #we take the df for the end residuals graph
  363. df_test = df_test.reset_index(drop=True)
  366. #%%Optional - grab demographics
  368. for dataset in datasets:
  369. df_clin = pd.read_csv(os.path.join(proc_dir, 'transfer_' + dataset, 'dem_patients.csv'))
  370. df_ctrl = pd.read_csv(os.path.join(proc_dir, 'transfer_' + dataset, 'dem_controls.csv'))
  372. print(f"{dataset} SCA patients")
  373. print(f"subjects total {len(df_clin.index.tolist())}")
  374. print(f"unique scanners {len(np.unique(df_clin["site"]))}")
  375. print(f"female {100*df_clin["sex"].isin(['F']).sum()/len(df_clin["sex"]):.1f} percent")
  376. print(f"male {100*df_clin["sex"].isin(['M']).sum()/len(df_clin["sex"]):.1f} percent")
  377. print(f"age mean {np.mean(df_clin["age"]):.1f}")
  378. print(f"age sd {np.std(df_clin["age"]):.1f}")
  379. print(f"SARA mean {np.mean(df_clin["SARA_score"]):.1f}")
  380. print(f"SARA sd {np.std(df_clin["SARA_score"]):.1f}")
  383. print(f"{dataset} controls")
  384. print(f"subjects total {len(df_ctrl.index.tolist())}")
  385. print(f"unique scanners {len(np.unique(df_ctrl["site"]))}")
  386. print(f"female {100*df_ctrl["sex"].isin(['F']).sum()/len(df_ctrl["sex"]):.1f} percent")
  387. print(f"male {100*df_ctrl["sex"].isin(['M']).sum()/len(df_ctrl["sex"]):.1f} percent")
  388. print(f"age mean {np.mean(df_ctrl["age"]):.1f}")
  389. print(f"age sd {np.std(df_ctrl["age"]):.1f}")
  391. #%% Configure prediction pipeline
  393. scaler = StandardScaler(with_mean = True, with_std = True)
  394. splitter = LeaveOneOut()
  395. inner_splitter = LeaveOneOut()
  397. #transformer to safely manage age residualization within pipeline
  398. class AgeResidualizer(BaseEstimator, TransformerMixin):
  399. def __init__(self, center_age=True):
  400. self.center_age = center_age
  402. def fit(self, X, y=None):
  404. age = X.iloc[:, -1].to_numpy().reshape(-1, 1)# last column of X is assumed to be age
  405. voxels = X.iloc[:, :-1].to_numpy()
  407. if self.center_age:
  408. self.age_mean_ = age.mean()
  409. age = age - self.age_mean_
  410. else:
  411. self.age_mean_ = 0.0
  413. Z = np.hstack([np.ones((age.shape[0], 1)), age])
  414. self.B_ = np.linalg.solve(Z.T @ Z, Z.T @ voxels)
  415. return self
  417. def transform(self, X):
  418. age = X.iloc[:, -1].to_numpy().reshape(-1, 1)
  419. voxels = X.iloc[:, :-1].to_numpy()
  421. age = age - self.age_mean_
  422. Z = np.hstack([np.ones((age.shape[0], 1)), age])
  423. voxels_res = voxels - Z @ self.B_
  425. return voxels_res #important: here we return only voxel columns and not the age one
  427. #we have to account for bias in KRR by adding an intercept column
  428. def add_intercept(X):
  429. intercept = np.ones((X.shape[0], 1))
  430. return np.hstack((X, intercept))
  432. add_intercept_transformer = FunctionTransformer(add_intercept)
  433. clf = KernelRidge()
  435. if options["Xstandard"] :
  436. pipeline = Pipeline(steps =[
  437. ('resid', AgeResidualizer()), #this pipeline expects age as last column of X_train and X_test and will remove it after
  438. ('scaler', scaler),
  439. ('intercept', add_intercept_transformer),
  440. ('clf', clf)])
  441. elif not options["Xstandard"] :
  442. pipeline = Pipeline(steps =[
  443. ('intercept', add_intercept_transformer),
  444. ('clf', clf)])
  446. #%% Run prediction
  448. fitted_models = []
  449. ys_pred = []
  450. ys_test = []
  451. ys_train = []
  452. Xs_train_final = []
  453. Xs_test_final = []
  454. feature_coefs = []
  455. eval_scores = pd.DataFrame(columns = ['RMSE'])
  457. for train_idx, test_idx in splitter.split(X_all):
  459. X_train, X_test = X_all.iloc[train_idx], X_all.iloc[test_idx]
  460. y_train, y_test = y_all[train_idx], y_all[test_idx]
  461. ys_test.append(y_test)
  462. ys_train.append(y_train)
  464. #fit
  465. pipeline.fit(X_train, y_train)
  466. best_pipeline = pipeline
  467. krr_model = best_pipeline.named_steps['clf']
  469. #extract outputs
  470. X_train_final = best_pipeline[:-1].transform(X_train) #extract scaled and possibly selected features
  471. X_test_final = best_pipeline[:-1].transform(X_test)
  472. Xs_train_final.append(X_train_final)
  473. Xs_test_final.append(X_test_final)
  475. coefs = krr_model.dual_coef_
  476. coefs = np.dot(coefs, X_train_final) #this only works for linear kernel of krr
  477. feature_coefs.append(coefs)
  479. y_pred = best_pipeline.predict(X_test)
  480. fitted_models.append(best_pipeline)
  481. ys_pred.append(y_pred)
  483. rmse = root_mean_squared_error(y_true = y_test, y_pred = y_pred)
  484. eval_scores.loc[len(eval_scores)] = [rmse]
  487. #%% Save outputs
  489. #keep syntax straight
  490. y_test = np.concatenate(ys_test)
  491. y_pred = np.concatenate(ys_pred)
  492. X_train_final = np.concatenate(Xs_train_final)
  493. X_test_final = np.concatenate(Xs_test_final)
  495. eval_scores['r2_score'] = [r2_score( y_true = y_test, y_pred = y_pred)] + [np.nan] * (len(eval_scores)-1)
  497. correlation_coef, p_value =pearsonr(y_test, y_pred)
  498. eval_scores['Pearson r'] = [correlation_coef] + [np.nan] * (len(eval_scores)-1)
  499. eval_scores['Pearson p-value'] = ["{:.2e}".format(p_value)] + [np.nan] * (len(eval_scores)-1)
  500. eval_scores.to_csv(os.path.join(final_dir,"_".join(['eval_scores']+[key for key, value in featureset.items() if value] +
  501. [key for key, value in train_and_test.items() if value] +
  502. [key for key, value in options.items() if value] )+ extension +'.csv'),
  503. index = False, na_rep = '')
  504. print(eval_scores.iloc[0,:])
  505. with open(os.path.join(final_dir,'fitted_models.pkl'), 'wb') as file:
  506. pickle.dump(fitted_models, file)
  507. with open(os.path.join(final_dir,'ys_pred.pkl'), 'wb') as file:
  508. pickle.dump(y_pred, file)
  509. with open(os.path.join(final_dir,'ys_test.pkl'), 'wb') as file:
  510. pickle.dump(y_test, file)
  511. with open(os.path.join(final_dir,'Xs_train.pkl'), 'wb') as file:
  512. pickle.dump(X_train_final, file)
  513. with open(os.path.join(final_dir,'Xs_test.pkl'), 'wb') as file:
  514. pickle.dump(X_test_final, file)
  515. with open(os.path.join(final_dir,'feature_coefs.pkl'), 'wb') as file:
  516. pickle.dump(feature_coefs, file)
  518. #map selected voxelwise feature mean coefs across folds back onto a nifti
  519. feature_coefs = [array[:-1] for array in feature_coefs] #remove last coefficient due to the added intercept (column of ones)
  520. average_coefs = np.mean(feature_coefs, axis = 0)
  521. ptksave(average_coefs, os.path.join(final_dir,'Signed_feature_coefs.nii.gz'), example=ex_nii, mask=mask_nii)
  523. #Plot residuals to check systematic over- or under-fitting
  524. diff_y = y_pred - y_test
  526. fig = plt.figure(figsize=(8, 6))
  527. y_name = ''.join(key for key, value in target.items() if value) #grab outcome name
  528. plt.title(f'{y_name}' + '(pred) - ' +f'{y_name}'+'(true)', fontsize=16)
  529. plt.plot(diff_y, marker='o', linestyle='')
  530. plt.gca().axes.xaxis.set_ticks([])
  531. plt.ylim(-20, 20)
  532. plt.axhline(y=0, color='gray', linestyle='--', linewidth=1)
  533. texts = []
  534. buffer = 0.7
  535. if target["Slope_BL_Y1_SARA"] : #adapt plotlims
  536. plt.ylim(-0.025, 0.025)
  537. buffer = 0.0008
  538. for i in range(len(diff_y)):
  539. text = plt.text(i, diff_y[i]+buffer,f'{df_test["sub_id"][i]}', ha='center', va='baseline',fontsize=8)
  540. plt.show()
  541. fig.savefig(os.path.join(final_dir,'y_pred-y_true.png'), dpi = 300, bbox_inches = 'tight')

06_KernelRidgeRegression_SCA_patients.py at commit ecb46eb, no license · at the source

Overview

Authors: Alice V. Chavanne1,2, Yaping Wang1,3,4, Augustijn A. A. de Boer1,2, Bing Xu5,6, Teije H. van Prooije1,7, Kirsten C. J. Kapteijns1,7, Colette Reniers1,7, Carlos R. Hernandez-Castillo8, Juan Fernandez-Ruiz9, Bart P. van de Warrenburg1,7, Jörn Diedrichsen3,4,10, Ryan L. Muetzel5,11, André F. Marquand1,2
ORCID iDs: Alice V. Chavanne, Bing Xu, Teije H. van Prooije, Kirsten C. J. Kapteijns, Colette Reniers, Carlos R. Hernandez-Castillo, Juan Fernandez-Ruiz, Bart P. van de Warrenburg, Jörn Diedrichsen, Ryan L. Muetzel, André F. Marquand
  1. Donders Institute for Brain, Cognition, and Behavior, Radboud University, Nijmegen, the Netherlands
  2. Department of Cognitive Neuroscience, Radboud University Medical Center, Nijmegen, the Netherlands
  3. Western Institute for Neuroscience, Western University, London, Ontario, Canada
  4. Department of Computer Science, Western University, London, Ontario, Canada
  5. Department of Child and Adolescent Psychology and Psychiatry, Erasmus MC, University Medical Center Rotterdam-Sophia Children ‘s Hospital, Rotterdam, The Netherlands
  6. The Generation R Study Group, Erasmus MC, University Medical Center Rotterdam, Rotterdam, The Netherlands
  7. Department of Neurology, Radboud University Medical Center, Nijmegen, the Netherlands
  8. Faculty of Computer Science, Dalhousie University, Halifax, Canada
  9. Departamento de Fisiología, Facultad de Medicina, Universidad Nacional Autónoma de México, Ciudad de Mexico, México
  10. Department of Statistical and Actuarial Sciences, Western University, London, Ontario, Canada
  11. Department of Radiology and Nuclear Medicine, Erasmus MC University Medical Center Rotterdam, Rotterdam, The Netherlands
Institutions: Donders Institute for Brain, Cognition and Behaviour (Netherlands); Radboud University Nijmegen (Netherlands); Radboud University Medical Center (Netherlands); Western University (Canada); Erasmus MC (Netherlands); Dalhousie University (Canada); Universidad Nacional Autónoma de México (Mexico)
Journal: bioRxiv (preprint)
Dates: published online 31 March 2026
Type: Preprint
License: CC BY
Identifiers: DOI 10.64898/2026.03.27.714738 · OpenAlex W7147445871
Open access: green, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism), other condition (population), computational (subfield)
Methods: Connectivity, Statistics, Machine learning, Preprocessing, fMRI & imaging, Physiology & signal measures
Topic: Genetic Neurodegenerative Diseases (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: European Commission (101057529); ZonMw (40‐44600‐98‐606); Raynor Cerebellum Project
Citations: not cited yet (Europe PMC); 64 references in the paper

Abstract

Brain disorders are often characterized by biological heterogeneity that is poorly captured by group-average analyses. Normative modeling has emerged as a promising tool to parse out such heterogeneity, yet existing lifespan reference models rely on coarse parcellations, which may obscure individual variability. Using an aggregated reference sample (n=58,597 scans from n=51,107 participants), we provide openly available normative models of brain morphometry at the voxel level across the lifespan, and we illustrate their potential utility with two complementary applications. First, we investigated long-term brain development after preterm birth across two independent cohorts (n=284; n=304) and found individualized, replicable and persistent brain alterations. Second, we extracted high-resolution patient-level morphometric deviations in two samples with rare, genetic neurodegenerative disorders (spinocerebellar ataxia type 1 and 3; n=29, n=15), which showed marked heterogeneity. Together, our findings highlight that voxelwise normative modeling can detect clinically relevant, individualized deviations from a reference model with high spatial precision.

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 12 matches between paragraphs and lines of code.

RaynorCerebellumCharts/NormativeModels

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: ecb46eb49e81de3679add0a77eaae910344502a9, 4 September 2026
Languages: Python (6), Jupyter (1)
Size: 18 files, 7 scripts
Software Heritage: not archived
Found in: the text, “Methods”
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (7 files), pandas (6 files), Matplotlib (5 files), NiBabel (5 files), scikit-learn (4 files), Nilearn (3 files), SciPy (3 files), seaborn (3 files), statsmodels (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
8 files
At the source: github.com/RaynorCerebellumCharts/NormativeModels

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

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, 30 September 2026: the first record

Recorded: type, journal, dates, 13 authors, 3 funders, 63 references.

Cite

This paper

Chavanne, A. V., Wang, Y., de Boer, A. A. A., Xu, B., van Prooije, T. H., Kapteijns, K. C. J., Reniers, C., Hernandez-Castillo, C. R., Fernandez-Ruiz, J., van de Warrenburg, B. P., Diedrichsen, J., Muetzel, R. L., & Marquand, A. F. (2026). Dissecting heterogeneous brain development and aging using voxelwise normative models. bioRxiv (preprint). https://doi.org/10.64898/2026.03.27.714738

BibTeX

@article{chavanne2026dissecting,
author = {Chavanne, Alice V. and Wang, Yaping and de Boer, Augustijn A. A. and Xu, Bing and van Prooije, Teije H. and Kapteijns, Kirsten C. J. and Reniers, Colette and Hernandez-Castillo, Carlos R. and Fernandez-Ruiz, Juan and van de Warrenburg, Bart P. and Diedrichsen, Jörn and Muetzel, Ryan L. and Marquand, André F.},
title = {{Dissecting heterogeneous brain development and aging using voxelwise normative models}},
journal = {bioRxiv (preprint)},
year = {2026},
month = mar,
publisher = {bioRxiv},
issn = {2692-8205},
doi = {10.64898/2026.03.27.714738},
url = {https://doi.org/10.64898/2026.03.27.714738}
}

RIS

TY - JOUR
AU - Chavanne, Alice V.
AU - Wang, Yaping
AU - de Boer, Augustijn A. A.
AU - Xu, Bing
AU - van Prooije, Teije H.
AU - Kapteijns, Kirsten C. J.
AU - Reniers, Colette
AU - Hernandez-Castillo, Carlos R.
AU - Fernandez-Ruiz, Juan
AU - van de Warrenburg, Bart P.
AU - Diedrichsen, Jörn
AU - Muetzel, Ryan L.
AU - Marquand, André F.
TI - Dissecting heterogeneous brain development and aging using voxelwise normative models
T2 - bioRxiv (preprint)
J2 - bioRxiv
PY - 2026
DA - 2026/03/31
SN - 2692-8205
PB - bioRxiv
DO - 10.64898/2026.03.27.714738
UR - https://doi.org/10.64898/2026.03.27.714738
ER -

CSL-JSON

{
"id": "10.64898/2026.03.27.714738",
"type": "article",
"title": "Dissecting heterogeneous brain development and aging using voxelwise normative models",
"container-title": "bioRxiv (preprint)",
"author": [
{
"family": "Chavanne",
"given": "Alice V."
},
{
"family": "Wang",
"given": "Yaping"
},
{
"family": "de Boer",
"given": "Augustijn A. A."
},
{
"family": "Xu",
"given": "Bing"
},
{
"family": "van Prooije",
"given": "Teije H."
},
{
"family": "Kapteijns",
"given": "Kirsten C. J."
},
{
"family": "Reniers",
"given": "Colette"
},
{
"family": "Hernandez-Castillo",
"given": "Carlos R."
},
{
"family": "Fernandez-Ruiz",
"given": "Juan"
},
{
"family": "van de Warrenburg",
"given": "Bart P."
},
{
"family": "Diedrichsen",
"given": "Jörn"
},
{
"family": "Muetzel",
"given": "Ryan L."
},
{
"family": "Marquand",
"given": "André F."
}
],
"container-title-short": "bioRxiv",
"DOI": "10.64898/2026.03.27.714738",
"ISSN": "2692-8205",
"publisher": "bioRxiv",
"URL": "https://doi.org/10.64898/2026.03.27.714738",
"issued": {
"date-parts": [
[
2026,
3,
31
]
]
}
}

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-024-46398-2 [code]
Population-wide cerebellar growth models of children and adolescents.
Journal: Nature communications
In common: NiBabel, seaborn, scikit-learn, 4 other tools, 7 references, 3 authors
[2] doi:10.1162/imag.a.1269 [code]
From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, statsmodels, NiBabel, 6 other tools, computational, 15 references
[3] doi:10.64898/2026.02.17.706268 [code]
Protocol Update: The Normative Modelling Paradigm for Computational Psychiatry
Journal: bioRxiv (preprint)
In common: statsmodels, NiBabel, seaborn, 5 other tools, computational, 14 references
[4] doi:10.1186/s13293-026-00846-4 [code]
Predicting developmental norms from baseline cortical thickness in longitudinal studies.
Journal: Biology of sex differences
In common: Nilearn, statsmodels, NiBabel, 6 other tools, 14 references
[5] doi:10.1038/s41467-026-72875-x [code]
Lifespan normative modeling of brain microstructure.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, 13 references, author Andre Marquand
[6] doi:10.1038/s42003-026-09825-2 [code]
Normative modeling of MEG brain oscillations across the human lifespan.
Journal: Communications biology
In common: Nilearn, statsmodels, NiBabel, 6 other tools, computational, 9 references, author Andre Marquand
[7] doi:10.1038/s42003-024-06573-z [code]
Dissecting task-based fMRI activity using normative modelling: an application to the Emotional Face Matching Task.
Journal: Communications biology
In common: statsmodels, NiBabel, seaborn, 4 other tools, 11 references, author Andre Marquand
[8] doi:10.1038/s41598-025-20484-x [code]
Understanding individual neurodegenerative progression in Parkinson's disease through normative modelling.
Journal: Scientific reports
In common: NiBabel, seaborn, scikit-learn, 4 other tools, computational, 12 references
[9] doi:10.1038/s41467-025-63967-1 [code]
Heterogeneous, temporally consistent, and plastic brain development after preterm birth.
Journal: Nature communications
In common: statsmodels, NiBabel, seaborn, 5 other tools, other condition, 11 references
[10] doi:10.1371/journal.pone.0278776 [code]
Closing the life-cycle of normative modeling using federated hierarchical Bayesian regression.
Journal: PloS one
In common: NiBabel, seaborn, scikit-learn, 4 other tools, computational, 14 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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.64898/2026.03.27.714738.svg)](https://openscicode.org/badge/doi/10.64898/2026.03.27.714738)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.64898/2026.03.27.714738.svg
   :target: https://openscicode.org/badge/doi/10.64898/2026.03.27.714738
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.64898/2026.03.27.714738"><img src="https://openscicode.org/badge/doi/10.64898/2026.03.27.714738.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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