OSCR

Exploring the link between body physiology and cognition: the role of the brain and aging.

Code ↔ Paper

20 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 20 matches
  1. [1] § Results › Predicting cognition from body phenotypes ↔ 9.Results/1Figures.ipynb, lines 1570–1656 · score 0.98 · peak expiratory flow, forced vital capacity, forced expiratory volume, muscle fat infiltration, head bone area, femur Ward
  2. [2] § Results › Predicting cognition from body phenotypes ↔ 9.Results/1Figures.ipynb, lines 1570–1656 · score 0.97 · systolic blood pressure, glycated haemoglobin, muscle fat infiltration, visceral adipose tissue, HbA1c, muscle ratio
  3. [3] § Methods › General cognition factor ↔ 1.Cognition/1.GetCognitiveScores.ipynb, lines 35–52 · score 0.95 · Symbol Digit Substitution, fluid intelligence, correctly identify matches, incorrect matches, complete numeric, Pairs Matching
  4. [4] § Results › Cognition–age relationship overlapping with body physiology and brain markers ↔ 6.Bootstrap/1Bootstrapping.ipynb, lines 467–611 · score 0.93 · Schaefer7n200p, ex vivo, aparc.a2009s, subcortical volumetric, FSL FAST, Schaefer7n500p
  5. [5] § Results › Cognition–age relationship overlapping with body physiology and brain markers ↔ 8.Commonality_mediation_analysis/1CommonalityAnalysis.ipynb, lines 265–412 · score 0.93 · Schaefer7n200p, ex vivo, aparc.a2009s, subcortical volumetric, FSL FAST, Schaefer7n500p
  6. [6] § Methods › Structural MRI (sMRI) ↔ 6.Bootstrap/1Bootstrapping.ipynb, lines 467–611 · score 0.91 · Desikan Killiany Tourville, subcortical volumetric subsegmentation, ex vivo, a2009s, FSL FAST, grey matter volumes
  7. [7] § Methods › Structural MRI (sMRI) ↔ 8.Commonality_mediation_analysis/1CommonalityAnalysis.ipynb, lines 265–412 · score 0.91 · Desikan Killiany Tourville, subcortical volumetric subsegmentation, ex vivo, a2009s, FSL FAST, grey matter volumes
  8. [8] § Methods › Diffusion-weighted MRI (dwMRI) ↔ 3_Get_brain_phenotypes/1GetMRIData_dwMRI.ipynb, lines 773–943 · score 0.91 · free water volume, diffusion tensor mode, orientation dispersion, probabilistic tractography, white matter tracts, volume fraction
  9. [9] § Methods › Body phenotypes ↔ 9.Results/2ResultSummary.ipynb, lines 74–108 · score 0.89 · pulse wave, carotid ultrasound, arterial stiffness, abdominal organs, heart MRI, electrocardiogram
  10. [10] § Methods › Body phenotypes ↔ 2.Get_body_phenotypes/1GetBodyPhenotypes.ipynb, lines 1522–1598 · score 0.89 · pulse wave, left ventricular, carotid ultrasound, arterial stiffness, heart MRI, fitness
  11. [11] § Results › Predicting cognition from body phenotypes ↔ 2.Get_body_phenotypes/1GetBodyPhenotypes.ipynb, lines 1437–1485 · score 0.87 · peak expiratory flow, forced vital capacity, forced expiratory volume, blood pressure, FEV1, FVC
  12. [12] § Results › Predicting cognition from body phenotypes ↔ 2.Get_body_phenotypes/1GetBodyPhenotypes.ipynb, lines 42–103 · score 0.84 · visceral adipose tissue, muscle fat infiltration, muscle ratio, fat fraction, thighs, pancreas
  13. [13] § Methods › Machine learning ↔ 5.Get_composite_markers_ML_2level_stacking/1-stack-body-brain.ipynb, lines 168–312 · score 0.69 · GridSearchCV, Random Forest, body brain, XGBoost, trained, stacking
  14. [14] § Results › Cognition–body relationship overlapping with brain markers ↔ 9.Results/1Figures.ipynb, lines 2673–2779 · score 0.69 · Venn diagrams, Stacked bar, Unique variance, common variance, individual body, composite
  15. [15] § Results ↔ 5.Get_composite_markers_ML_2level_stacking/1-stack-body-brain.ipynb, lines 168–312 · score 0.63 · Random Forest, absolute error, XGBoost, CV, algorithm, MAE
  16. [16] § Methods › Resting-state functional MRI (rsMRI) ↔ 3_Get_brain_phenotypes/2GetMRIData_rsMRI_Transpose_GetCorrelationNilearn.py, lines 17–119 · score 0.61 · ConnectivityMeasure, partial correlation matrices, Nilearn, atlas, brain phenotypes, MRI
  17. [17] § Methods › Resting-state functional MRI (rsMRI) ↔ 3_Get_brain_phenotypes/1GetMRIData_rsMRI.ipynb, lines 465–608 · score 0.61 · ConnectivityMeasure, partial correlation matrices, Nilearn, rsMRI, atlas, brain phenotypes
  18. [18] § Methods › Machine learning ↔ 5.Get_composite_markers_ML_2level_stacking/1-stack-body-brain.ipynb, lines 19–166 · score 0.60 · absolute error, squared error, configuration, MSE, MAE, outer
  19. [19] § Results › Predicting cognition from brain phenotypes ↔ 3_Get_brain_phenotypes/5GetMRIData_sMRI.ipynb, lines 251–424 · score 0.58 · choroid plexus, lateral ventricle, hypointensities, cortex, posterior, brain phenotypes
  20. [20] § Methods › Neuroimaging ↔ 3_Get_brain_phenotypes/1GetMRIData_dwMRI.ipynb, lines 773–943 · score 0.56 · fiber orientation, water diffusion, dwMRI, axons, weighted, density

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 · 3,791 lines · 137 KB · MIT · 3 matches

  1. # %%
  2. import csv
  3. import os
  4. import re
  5. import random
  6. import pickle
  7. import openpyxl
  8. import warnings
  9. import pandas as pd
  10. import scipy
  11. import numpy as np
  12. import matplotlib as mpl
  13. from matplotlib_venn import venn2
  14. from matplotlib_venn import venn3
  15. import matplotlib.gridspec as gridspec
  16. from matplotlib.gridspec import GridSpec
  17. from matplotlib import ticker
  18. import matplotlib.ticker as mticker
  19. from matplotlib.legend_handler import HandlerPatch
  20. import matplotlib.cm as cm
  21. import matplotlib.patches as mpatches
  22. import matplotlib.colors as colors
  23. from scipy import stats
  24. from scipy.stats import gaussian_kde
  25. import statsmodels.api as sm
  26. import matplotlib.pyplot as plt
  27. from matplotlib.ticker import FormatStrFormatter, MultipleLocator, AutoMinorLocator, MaxNLocator, FixedLocator
  28. import seaborn as sns
  29. from scipy.stats import pearsonr
  30. from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
  31. from matplotlib.collections import PolyCollection
  32. import matplotlib.patches as patches
  33. from tqdm import tqdm
  34. from typing import List, Dict, Optional
  36. # %%
  37. # Define modality renaming dictionary (modality_map)
  38. modality_map = {
  39. 'hearing': 'Hearing',
  40. 'immune': 'Immune',
  41. 'renalhepatic': 'Renal & Hepatic',
  42. 'metabolic': 'Metabolic',
  43. 'cardiopulmonary': 'Cardiopulmonary',
  44. 'musculoskeletal': 'Musculoskeletal',
  45. 'bone_densitometry': 'Bone Densitometry of Heel',
  46. 'pwa': 'Pulse Wave Analysis',
  47. 'heart_mri': 'Heart MRI',
  48. 'carotid_ultrasound': 'Carotid Ultrasound',
  49. 'arterial_stiffness': 'Arterial Stiffness',
  50. 'ecg_rest': 'ECG at Rest',
  51. 'body_composition_by_impedance': 'Body Composition by Impedance',
  52. 'body_composition_dxa': 'Body Composition by DXA',
  53. 'bone_dxa': 'Bone Size, Mineral and Density by DXA',
  54. 'kidneys_mri': 'Kidney MRI',
  55. 'liver_mri': 'Liver MRI',
  56. 'abdominal_composition_mri_18_vars': 'Abdominal Composition by MRI',
  57. 'abdominal_organ_composition_mri_13_vars': 'Abdominal Organ Composition by MRI',
  58. 'struct_fast' : 'Regional grey matter volumes (FSL FAST)',
  59. 'struct_sub_first': 'Subcortical volumes (FSL FIRST)',
  61. 'struct_fs_aseg_mean_intensity' : 'ASEG Mean Intensity',
  62. 'struct_fs_aseg_volume' : 'ASEG Volume',
  65. 'struct_ba_exvivo_area' : 'BA ex-vivo Area',
  66. 'struct_ba_exvivo_mean_thickness' : 'BA ex-vivo Mean Thickness',
  67. 'struct_ba_exvivo_volume' : 'BA ex-vivo Volume',
  69. 'struct_a2009s_area' : 'a2009s Area',
  70. 'struct_a2009s_mean_thickness' : 'a2009s Mean Thickness',
  71. 'struct_a2009s_volume' : 'a2009s Volume',
  74. 'struct_dkt_area' : 'Desikan-Killiany-Tourville Area',
  75. 'struct_dkt_mean_thickness' : 'Desikan-Killiany-Tourville Mean Thickness',
  76. 'struct_dkt_volume' : 'Desikan-Killiany-Tourville Volume',
  79. 'struct_desikan_gw' : 'Desikan Grey/White Matter Contrast',
  80. 'struct_desikan_pial' : 'Desikan Pial',
  82. 'struct_desikan_white_area' : 'Desikan White Matter Area',
  83. 'struct_desikan_white_mean_thickness' : 'Desikan White Matter Mean Thickness',
  84. 'struct_desikan_white_volume' : 'Desikan White Matter Volume',
  85. "struct_subsegmentation":'Subcortical Volumetric Subsegmentation',
  87. 'add_t1' : 'Whole-Brain T1w',
  88. 'add_t2' : 'Whole-Brain T2w',
  89. "dwi_FA_tbss": "FA TBSS",
  90. "dwi_FA_prob": "FA Probabilistic",
  91. "dwi_MD_tbss": "MD TBSS",
  92. "dwi_MD_prob": "MD Probabilistic",
  93. "dwi_L1_tbss": "L1 TBSS",
  94. "dwi_L1_prob": "L1 Probabilistic",
  95. "dwi_L2_tbss": "L2 TBSS",
  96. "dwi_L2_prob": "L2 Probabilistic",
  97. "dwi_L3_tbss": "L3 TBSS",
  98. "dwi_L3_prob": "L3 Probabilistic",
  99. "dwi_MO_tbss": "MO TBSS",
  100. "dwi_MO_prob": "MO Probabilistic",
  101. "dwi_OD_tbss": "OD TBSS",
  102. "dwi_OD_prob": "OD Probabilistic",
  103. "dwi_ICVF_tbss": "ICVF TBSS",
  104. "dwi_ICVF_prob": "ICVF Probabilistic",
  105. "dwi_ISOVF_tbss": "ISOVF TBSS",
  106. "dwi_ISOVF_prob": 'ISOVF Probabilistic',
  107. "amplitudes_21": " 21 IC Amplitudes",
  108. "amplitudes_55": "55 IC Amplitudes",
  109. "full_correlation_21": "21 IC Full Correlation",
  110. "full_correlation_55": "55 IC Full Correlation",
  111. "partial_correlation_21": " 21 IC Partial Correlation",
  112. "partial_correlation_55": " 55 IC Partial Correlation",
  113. # aparc Tian S1 (I)
  114. 'aparc_Tian_S1_FA_i2': 'aparc-I FA',
  115. 'aparc_Tian_S1_Length_i2': 'aparc-I Length',
  116. 'aparc_Tian_S1_SIFT2_FBC_i2': 'aparc-I SIFT2 FBC',
  117. 'aparc_Tian_S1_Streamline_Count_i2': 'aparc-I Streamline Count',
  119. # aparc a2009s Tian S1 (I)
  120. 'aparc_a2009s_Tian_S1_FA_i2': 'aparc.a2009s-I FA',
  121. 'aparc_a2009s_Tian_S1_Length_i2': 'aparc.a2009s-I Length',
  122. 'aparc_a2009s_Tian_S1_SIFT2_FBC_i2': 'aparc.a2009s-I SIFT2 FBC',
  123. 'aparc_a2009s_Tian_S1_Streamline_Count_i2': 'aparc.a2009s-I Streamline Count',
  125. # Glasser Tian S1 (I)
  126. 'Glasser_Tian_S1_FA_i2': 'Glasser-I FA',
  127. 'Glasser_Tian_S1_Length_i2': 'Glasser-I Length',
  128. 'Glasser_Tian_S1_SIFT2_FBC_i2': 'Glasser-I SIFT2 FBC',
  129. 'Glasser_Tian_S1_Streamline_Count_i2': 'Glasser-I Streamline Count',
  131. # Glasser Tian S4 (IV)
  132. 'Glasser_Tian_S4_FA_i2': 'Glasser-IV FA',
  133. 'Glasser_Tian_S4_Length_i2': 'Glasser-IV Length',
  134. 'Glasser_Tian_S4_SIFT2_FBC_i2': 'Glasser-IV SIFT2 FBC',
  135. 'Glasser_Tian_S4_Streamline_Count_i2': 'Glasser-IV Streamline Count',
  137. # Schaefer7n1000p Tian S4 (IV) (in reality: Schaefer7n200p Tian S1)
  138. 'Schaefer7n1000p_Tian_S4_FA_i2': 'Schaefer7n200p-I FA', #'Schaefer7n1000p-IV FA',
  139. 'Schaefer7n1000p_Tian_S4_Length_i2': 'Schaefer7n200p-I Length',#'Schaefer7n1000p-IV Length',
  140. 'Schaefer7n1000p_Tian_S4_SIFT2_FBC_i2': 'Schaefer7n200p-I SIFT2 FBC',#'Schaefer7n1000p-IV SIFT2 FBC',
  141. 'Schaefer7n1000p_Tian_S4_Streamline_Count_i2': 'Schaefer7n200p-I Streamline Count', #'Schaefer7n1000p-IV Streamline Count'
  143. # Schaefer7n200p Tian S4 (IV) (in reality: Schaefer7n500p Tian S4)
  144. 'Schaefer7n200p_Tian_S1_FA_i2': 'Schaefer7n500p-IV FA',
  145. 'Schaefer7n200p_Tian_S1_Length_i2': 'Schaefer7n500p-IV Length',
  146. 'Schaefer7n200p_Tian_S1_SIFT2_FBC_i2': 'Schaefer7n500p-IV SIFT2 FBC',
  147. 'Schaefer7n200p_Tian_S1_Streamline_Count_i2': 'Schaefer7n500p-IV Streamline Count',
  149. # Schaefer7n500p Tian S4 (IV) (in reality: Schaefer7n1000p Tian S4)
  150. 'Schaefer7n500p_Tian_S4_FA_i2': 'Schaefer7n1000p-IV FA',
  151. 'Schaefer7n500p_Tian_S4_Length_i2': 'Schaefer7n1000p-IV Length',
  152. 'Schaefer7n500p_Tian_S4_SIFT2_FBC_i2': 'Schaefer7n1000p-IV SIFT2 FBC',
  153. 'Schaefer7n500p_Tian_S4_Streamline_Count_i2': 'Schaefer7n1000p-IV Streamline Count',
  155. # Resting state
  156. 'full_correlation_aparc_a2009s_Tian_S1' : 'aparc.a2009s-I Full Correlation',
  157. 'full_correlation_aparc_Tian_S1': 'aparc-I Full Correlation',
  158. 'full_correlation_Glasser_Tian_S1': 'Glasser-I Full Correlation',
  159. 'full_correlation_Glasser_Tian_S4': 'Glasser-IV Full Correlation',
  160. 'full_correlation_Schaefer7n200p_Tian_S1': 'Schaefer7n200p-I Full Correlation',
  161. 'full_correlation_Schaefer7n500p_Tian_S4': 'Schaefer7n500p-IV Full Correlation',
  162. 'partial_correlation_aparc_a2009s_Tian_S1': 'aparc.a2009s-I Partial Correlation',
  163. 'partial_correlation_aparc_Tian_S1': 'aparc-I Partial Correlation',
  164. 'partial_correlation_Glasser_Tian_S1': 'Glasser-I Partial Correlation',
  165. 'partial_correlation_Glasser_Tian_S4': 'Glasser-IV Partial Correlation',
  166. 'partial_correlation_Schaefer7n200p_Tian_S1': 'Schaefer7n200p-I Partial Correlation',
  167. 'partial_correlation_Schaefer7n500p_Tian_S4': 'Schaefer7n500p-IV Partial Correlation',
  170. 'allmri': '3 Brain MRI Modalities Stacked',
  171. 'dwi': 'Brain dwMRI Stacked',
  172. 'smri': 'Brain sMRI Stacked',
  173. 'rs': 'Brain rsMRI Stacked',
  174. 'body': 'Body Physiology Stacked',
  175. 'brain-plus-body': '3 Brain MRI Modalities & Body Stacked',
  176. 'brain-body': 'Brain & Body Stacked',
  177. 'body-only': 'Body Physiology and Composition Stacked'
  178. }
  180. # Define modality names for renaming
  181. modality_names = {
  182. 'hearing': 'Hearing',
  183. 'immune': 'Immune',
  184. 'renalhepatic': 'Renal & Hepatic',
  185. 'metabolic': 'Metabolic',
  186. 'cardiopulmonary': 'Cardiopulmonary',
  187. 'musculoskeletal': 'Musculoskeletal',
  188. 'bone_densitometry': 'Bone Densitometry of Heel',
  189. 'pwa': 'Pulse Wave Analysis',
  190. 'heart_mri': 'Heart MRI',
  191. 'carotid_ultrasound': 'Carotid Ultrasound',
  192. 'arterial_stiffness': 'Arterial Stiffness',
  193. 'ecg_rest': 'ECG at Rest',
  194. 'body_composition_by_impedance': 'Body Composition by Impedance',
  195. 'body_composition_dxa': 'Body Composition by DXA',
  196. 'bone_dxa': 'Bone Size, Mineral and Density by DXA',
  197. 'kidneys_mri': 'Kidney MRI',
  198. 'liver_mri': 'Liver MRI',
  199. 'abdominal_composition_mri_18_vars': 'Abdominal Composition by MRI',
  200. 'abdominal_organ_composition_mri_13_vars': 'Abdominal Organ Composition by MRI',
  201. 'struct_fast' : 'Regional grey matter volumes (FSL FAST)',
  202. 'struct_sub_first': 'Subcortical volumes (FSL FIRST)',
  204. 'struct_fs_aseg_mean_intensity' : 'ASEG Mean Intensity',
  205. 'struct_fs_aseg_volume' : 'ASEG Volume',
  208. 'struct_ba_exvivo_area' : 'BA ex-vivo Area',
  209. 'struct_ba_exvivo_mean_thickness' : 'BA ex-vivo Mean Thickness',
  210. 'struct_ba_exvivo_volume' : 'BA ex-vivo Volume',
  212. 'struct_a2009s_area' : 'a2009s Area',
  213. 'struct_a2009s_mean_thickness' : 'a2009s Mean Thickness',
  214. 'struct_a2009s_volume' : 'a2009s Volume',
  217. 'struct_dkt_area' : 'Desikan-Killiany-Tourville Area',
  218. 'struct_dkt_mean_thickness' : 'Desikan-Killiany-Tourville Mean Thickness',
  219. 'struct_dkt_volume' : 'Desikan-Killiany-Tourville Volume',
  222. 'struct_desikan_gw' : 'Desikan Grey/White Matter Contrast',
  223. 'struct_desikan_pial' : 'Desikan Pial',
  225. 'struct_desikan_white_area' : 'Desikan White Matter Area',
  226. 'struct_desikan_white_mean_thickness' : 'Desikan White Matter Mean Thickness',
  227. 'struct_desikan_white_volume' : 'Desikan White Matter Volume',
  228. "struct_subsegmentation":'Subcortical Volumetric Subsegmentation',
  230. 'add_t1' : 'Whole-Brain T1w',
  231. 'add_t2' : 'Whole-Brain T2w',
  232. "dwi_FA_tbss": "FA TBSS",
  233. "dwi_FA_prob": "FA Probabilistic",
  234. "dwi_MD_tbss": "MD TBSS",
  235. "dwi_MD_prob": "MD Probabilistic",
  236. "dwi_L1_tbss": "L1 TBSS",
  237. "dwi_L1_prob": "L1 Probabilistic",
  238. "dwi_L2_tbss": "L2 TBSS",
  239. "dwi_L2_prob": "L2 Probabilistic",
  240. "dwi_L3_tbss": "L3 TBSS",
  241. "dwi_L3_prob": "L3 Probabilistic",
  242. "dwi_MO_tbss": "MO TBSS",
  243. "dwi_MO_prob": "MO Probabilistic",
  244. "dwi_OD_tbss": "OD TBSS",
  245. "dwi_OD_prob": "OD Probabilistic",
  246. "dwi_ICVF_tbss": "ICVF TBSS",
  247. "dwi_ICVF_prob": "ICVF Probabilistic",
  248. "dwi_ISOVF_tbss": "ISOVF TBSS",
  249. "dwi_ISOVF_prob": 'ISOVF Probabilistic',
  250. "amplitudes_21": " 21 IC Amplitudes",
  251. "amplitudes_55": "55 IC Amplitudes",
  252. "full_correlation_21": "21 IC Full Correlation",
  253. "full_correlation_55": "55 IC Full Correlation",
  254. "partial_correlation_21": " 21 IC Partial Correlation",
  255. "partial_correlation_55": " 55 IC Partial Correlation",
  256. # aparc Tian S1 (I)
  257. 'aparc_Tian_S1_FA_i2': 'aparc-I FA',
  258. 'aparc_Tian_S1_Length_i2': 'aparc-I Length',
  259. 'aparc_Tian_S1_SIFT2_FBC_i2': 'aparc-I SIFT2 FBC',
  260. 'aparc_Tian_S1_Streamline_Count_i2': 'aparc-I Streamline Count',
  262. # aparc a2009s Tian S1 (I)
  263. 'aparc_a2009s_Tian_S1_FA_i2': 'aparc.a2009s-I FA',
  264. 'aparc_a2009s_Tian_S1_Length_i2': 'aparc.a2009s-I Length',
  265. 'aparc_a2009s_Tian_S1_SIFT2_FBC_i2': 'aparc.a2009s-I SIFT2 FBC',
  266. 'aparc_a2009s_Tian_S1_Streamline_Count_i2': 'aparc.a2009s-I Streamline Count',
  268. # Glasser Tian S1 (I)
  269. 'Glasser_Tian_S1_FA_i2': 'Glasser-I FA',
  270. 'Glasser_Tian_S1_Length_i2': 'Glasser-I Length',
  271. 'Glasser_Tian_S1_SIFT2_FBC_i2': 'Glasser-I SIFT2 FBC',
  272. 'Glasser_Tian_S1_Streamline_Count_i2': 'Glasser-I Streamline Count',
  274. # Glasser Tian S4 (IV)
  275. 'Glasser_Tian_S4_FA_i2': 'Glasser-IV FA',
  276. 'Glasser_Tian_S4_Length_i2': 'Glasser-IV Length',
  277. 'Glasser_Tian_S4_SIFT2_FBC_i2': 'Glasser-IV SIFT2 FBC',
  278. 'Glasser_Tian_S4_Streamline_Count_i2': 'Glasser-IV Streamline Count',
  280. # Schaefer7n1000p Tian S4 (IV) (in reality: Schaefer7n200p Tian S1)
  281. 'Schaefer7n1000p_Tian_S4_FA_i2': 'Schaefer7n200p-I FA', #'Schaefer7n1000p-IV FA',
  282. 'Schaefer7n1000p_Tian_S4_Length_i2': 'Schaefer7n200p-I Length',#'Schaefer7n1000p-IV Length',
  283. 'Schaefer7n1000p_Tian_S4_SIFT2_FBC_i2': 'Schaefer7n200p-I SIFT2 FBC',#'Schaefer7n1000p-IV SIFT2 FBC',
  284. 'Schaefer7n1000p_Tian_S4_Streamline_Count_i2': 'Schaefer7n200p-I Streamline Count', #'Schaefer7n1000p-IV Streamline Count'
  286. # Schaefer7n200p Tian S4 (IV) (in reality: Schaefer7n500p Tian S4)
  287. 'Schaefer7n200p_Tian_S1_FA_i2': 'Schaefer7n500p-IV FA',
  288. 'Schaefer7n200p_Tian_S1_Length_i2': 'Schaefer7n500p-IV Length',
  289. 'Schaefer7n200p_Tian_S1_SIFT2_FBC_i2': 'Schaefer7n500p-IV SIFT2 FBC',
  290. 'Schaefer7n200p_Tian_S1_Streamline_Count_i2': 'Schaefer7n500p-IV Streamline Count',
  292. # Schaefer7n500p Tian S4 (IV) (in reality: Schaefer7n1000p Tian S4)
  293. 'Schaefer7n500p_Tian_S4_FA_i2': 'Schaefer7n1000p-IV FA',
  294. 'Schaefer7n500p_Tian_S4_Length_i2': 'Schaefer7n1000p-IV Length',
  295. 'Schaefer7n500p_Tian_S4_SIFT2_FBC_i2': 'Schaefer7n1000p-IV SIFT2 FBC',
  296. 'Schaefer7n500p_Tian_S4_Streamline_Count_i2': 'Schaefer7n1000p-IV Streamline Count',
  298. # Resting state
  299. 'full_correlation_aparc_a2009s_Tian_S1' : 'aparc.a2009s-I Full Correlation',
  300. 'full_correlation_aparc_Tian_S1': 'aparc-I Full Correlation',
  301. 'full_correlation_Glasser_Tian_S1': 'Glasser-I Full Correlation',
  302. 'full_correlation_Glasser_Tian_S4': 'Glasser-IV Full Correlation',
  303. 'full_correlation_Schaefer7n200p_Tian_S1': 'Schaefer7n200p-I Full Correlation',
  304. 'full_correlation_Schaefer7n500p_Tian_S4': 'Schaefer7n500p-IV Full Correlation',
  305. 'partial_correlation_aparc_a2009s_Tian_S1': 'aparc.a2009s-I Partial Correlation',
  306. 'partial_correlation_aparc_Tian_S1': 'aparc-I Partial Correlation',
  307. 'partial_correlation_Glasser_Tian_S1': 'Glasser-I Partial Correlation',
  308. 'partial_correlation_Glasser_Tian_S4': 'Glasser-IV Partial Correlation',
  309. 'partial_correlation_Schaefer7n200p_Tian_S1': 'Schaefer7n200p-I Partial Correlation',
  310. 'partial_correlation_Schaefer7n500p_Tian_S4': 'Schaefer7n500p-IV Partial Correlation',
  312. 'lifestyle-envir': 'Lifestyle & Environment',
  314. 'allmri': '3 Brain MRI Modalities Stacked',
  315. 'dwi': 'Brain dwMRI Stacked',
  316. 'smri': 'Brain sMRI Stacked',
  317. 'rs': 'Brain rsMRI Stacked',
  318. 'body': 'Body Physiology Stacked',
  319. 'body-comp': 'Body Composition Stacked',
  320. #'cardiopulmonary': 'Cardiopulmonary stacked',
  321. 'renal-hepatic': 'Renal & Hepatic Stacked',
  322. 'lifestyle-envir': 'Lifestyle & Environment',
  323. 'brain-plus-body': '3 Brain MRI Modalities & Body Stacked',
  324. 'brain-body': 'Brain & Body Stacked',
  325. 'body-only': 'Body Physiology and Composition Stacked'
  326. }
  328. # %%
  329. # Define brain and body modalities
  330. modalities_smri = [
  331. 'struct_fast',
  332. 'struct_sub_first',
  333. 'struct_fs_aseg_mean_intensity',
  334. 'struct_fs_aseg_volume',
  335. 'struct_ba_exvivo_area',
  336. 'struct_ba_exvivo_mean_thickness',
  337. 'struct_ba_exvivo_volume',
  338. 'struct_a2009s_area',
  339. 'struct_a2009s_mean_thickness',
  340. 'struct_a2009s_volume',
  341. 'struct_dkt_area',
  342. 'struct_dkt_mean_thickness',
  343. 'struct_dkt_volume',
  344. 'struct_desikan_gw',
  345. 'struct_desikan_pial',
  346. 'struct_desikan_white_area',
  347. 'struct_desikan_white_mean_thickness',
  348. 'struct_desikan_white_volume',
  349. 'struct_subsegmentation',
  350. 'add_t1',
  351. 'add_t2']
  353. modalities_dwmri = [
  355. "dwi_FA_tbss", "dwi_FA_prob",
  356. "dwi_MD_tbss", "dwi_MD_prob",
  357. "dwi_L1_tbss", "dwi_L1_prob",
  358. "dwi_L2_tbss", "dwi_L2_prob",
  359. "dwi_L3_tbss", "dwi_L3_prob",
  360. "dwi_MO_tbss", "dwi_MO_prob",
  361. "dwi_OD_tbss", "dwi_OD_prob",
  362. "dwi_ICVF_tbss", "dwi_ICVF_prob",
  363. "dwi_ISOVF_tbss", "dwi_ISOVF_prob",
  365. 'aparc_Tian_S1_FA_i2',
  366. 'aparc_Tian_S1_Length_i2',
  367. 'aparc_Tian_S1_SIFT2_FBC_i2',
  368. 'aparc_Tian_S1_Streamline_Count_i2',
  370. 'aparc_a2009s_Tian_S1_FA_i2',
  371. 'aparc_a2009s_Tian_S1_Length_i2',
  372. 'aparc_a2009s_Tian_S1_SIFT2_FBC_i2',
  373. 'aparc_a2009s_Tian_S1_Streamline_Count_i2',
  375. 'Glasser_Tian_S1_FA_i2',
  376. 'Glasser_Tian_S1_Length_i2',
  377. 'Glasser_Tian_S1_SIFT2_FBC_i2',
  378. 'Glasser_Tian_S1_Streamline_Count_i2',
  380. 'Glasser_Tian_S4_FA_i2',
  381. 'Glasser_Tian_S4_Length_i2',
  382. 'Glasser_Tian_S4_SIFT2_FBC_i2',
  383. 'Glasser_Tian_S4_Streamline_Count_i2',
  385. 'Schaefer7n200p_Tian_S1_FA_i2',
  386. 'Schaefer7n200p_Tian_S1_Length_i2',
  387. 'Schaefer7n200p_Tian_S1_SIFT2_FBC_i2',
  388. 'Schaefer7n200p_Tian_S1_Streamline_Count_i2',
  390. 'Schaefer7n1000p_Tian_S4_FA_i2',
  391. 'Schaefer7n1000p_Tian_S4_Length_i2',
  392. 'Schaefer7n1000p_Tian_S4_SIFT2_FBC_i2',
  393. 'Schaefer7n1000p_Tian_S4_Streamline_Count_i2']
  395. modalities_rsmri = [
  397. "amplitudes_21",
  398. "full_correlation_21",
  399. "partial_correlation_21",
  400. "amplitudes_55",
  401. "full_correlation_55",
  402. "partial_correlation_55",
  403. 'full_correlation_aparc_a2009s_Tian_S1',
  404. 'full_correlation_aparc_Tian_S1',
  405. 'full_correlation_Glasser_Tian_S1',
  406. 'full_correlation_Glasser_Tian_S4',
  407. 'full_correlation_Schaefer7n200p_Tian_S1',
  408. 'full_correlation_Schaefer7n500p_Tian_S4',
  409. 'partial_correlation_aparc_a2009s_Tian_S1',
  410. 'partial_correlation_aparc_Tian_S1',
  411. 'partial_correlation_Glasser_Tian_S1',
  412. 'partial_correlation_Glasser_Tian_S4',
  413. 'partial_correlation_Schaefer7n200p_Tian_S1',
  414. 'partial_correlation_Schaefer7n500p_Tian_S4'
  415. ]
  417. modalities_body = [
  418. 'immune',
  419. 'renalhepatic',
  420. 'metabolic',
  421. 'cardiopulmonary',
  422. 'musculoskeletal',
  423. 'bone_densitometry',
  424. 'pwa',
  425. 'heart_mri',
  426. 'carotid_ultrasound',
  427. 'arterial_stiffness',
  428. 'ecg_rest',
  429. 'body_composition_by_impedance',
  430. 'body_composition_dxa',
  431. 'bone_dxa',
  432. 'kidneys_mri',
  433. 'liver_mri',
  434. 'abdominal_composition_mri_18_vars', #17 vars
  435. 'abdominal_organ_composition_mri_13_vars', #12 vars
  436. 'hearing'
  437. ]
  440. # %%
  441. # Define figures path
  442. fig_path = '/UK_BB/brainbody/figures'
  444. # %%
  445. # Define a function to build circle legend handles
  446. from matplotlib.legend_handler import HandlerPatch
  447. def make_circle_legend(legend, orig_handle, xdescent, ydescent, width, height, fontsize):
  448. # Create a circle for the legend
  449. return plt.Circle((width/2, height/2), min(width, height))
  451. # %%
  452. # Define a function to save figure in multiple formats
  453. def save_fig(formats, path):
  454. for fmt in formats:
  455. plt.savefig(
  456. f'{path}.{fmt}',
  457. bbox_inches="tight",
  458. pad_inches=1,
  459. transparent=False,
  460. facecolor="w",
  461. edgecolor='w',
  462. orientation='landscape',
  463. format=fmt,
  464. dpi=300
  465. )
  467. # %% [markdown]
  468. # # Figure 1
  470. # %% [markdown]
  471. # g-factor loading plot
  473. # %%
  474. # Loading plot
  475. path = "/UK_BB/brainbody/cognition/code/efa_loadings/efa_loadings_renamed.xlsx"
  476. df = pd.read_excel(path, sheet_name="Group_0")
  478. # Extract variable names
  479. variables = df["Variable"].astype(str).tolist()
  481. # Extract the four components
  482. components = ["PA1", "PA2", "PA3", "PA4"]
  483. title_map = {'PA1': 'Factor 1',
  484. 'PA2': 'Factor 2',
  485. 'PA3': 'Factor 3',
  486. 'PA4': 'Factor 4'}
  488. label_map = {
  489. 'Numeric memory: Maximum digits remembered correctly': 'Numeric memory: Max digits correct',
  490. 'Prospective memory: Initial answer': 'Prospective memory: Initial',
  491. 'Matrix pattern completion: Number of puzzles correct': 'Matrix pattern: Number correct',
  492. 'Fluid intelligence score': 'FIS',
  493. 'PAL: Number of correct pairs': 'PAL: Pairs correct',
  494. 'Tower rearranging: Number of puzzles correct': 'Tower rearranging: Puzzles correct',
  495. 'Picture vocabulary: Specific cognitive ability': 'Picture vocab.: Specific cog. ability',
  496. 'Symbol digit substitution: Proportion of correct matches': 'SDS: Proportion correct',
  497. '(log)TMT: Duration to complete numeric path': '(log)TMT: Duration numeric',
  498. '(log)TMT: Duration to complete alphabetic path': '(log)TMT: Duration alphabetic',
  499. '(log)Reaction time': '(log)RT',
  500. '(logx+1)Pairs matching: Incorrect matches': '(logx+1)Pairs matching: Incorrect'
  501. }
  502. renamed_variables = [label_map.get(v, v) for v in variables]
  504. # 2x2 subplots
  505. fig, axes = plt.subplots(2, 2, figsize=(6, 8), sharey=True)
  507. # Flatten axes so we can iterate cleanly
  508. axes = axes.flatten()
  510. for ax, comp in zip(axes, components):
  511. loadings = df[comp].astype(float).tolist()
  512. colors = ["#4A7169" if v >= 0 else "#A8554E" for v in loadings]
  514. ax.barh(variables, loadings, color=colors)
  515. ax.axvline(0, color="black", linewidth=1)
  517. all_loadings = pd.concat([df[c] for c in components])
  518. xmin = all_loadings.min()
  519. xmax = all_loadings.max()
  520. ax.set_xlim(xmin, xmax)
  521. ax.xaxis.set_major_locator(MultipleLocator(0.2))
  523. ax.set_title(title_map[comp], fontsize=22)
  524. ax.set_yticks(range(len(variables)))
  526. # Position tick labels
  527. if ax in axes:
  528. ax.yaxis.set_ticks_position('left')
  529. ax.tick_params(length=3)
  530. ax.tick_params(axis='y', length=2)
  532. # Set y-tick labels to show 1-12 from top to bottom
  533. ax.set_yticklabels(range(12, 0, -1), fontsize=17)
  535. ax.tick_params(axis='x', labelsize=18)
  536. ax.set_xlabel("Loading", fontsize=19)
  538. plt.tight_layout()
  539. plt.savefig(os.path.join(fig_path, 'fig1-prep/Fig1b_g_loadings.svg'),
  540. bbox_inches="tight",
  541. transparent=False,
  542. facecolor="w",
  543. edgecolor='w',
  544. orientation='landscape',
  545. format='svg',
  546. dpi=300)
  547. plt.show()
  549. # Print variable names in the order they appear in the plot (top to bottom)
  550. print("Variables in plot:")
  551. print("=" * 60)
  552. for i, var in enumerate(reversed(renamed_variables), 1):
  553. print(f"{13-i}. {var}")
  555. # %% [markdown]
  556. # ## Figure 2
  558. # %% [markdown]
  559. # Body performance
  561. # %%
  562. # KDE plot: Body composition by impedance vs composite body + Ridge plots for body performance
  563. # First row: KDE plot (1/3 height), Second row: violin plots (2/3 height)
  564. warnings.simplefilter(action='ignore', category=FutureWarning)
  566. fig = plt.figure(figsize=(14, 15))
  567. gs = fig.add_gridspec(2, 1, height_ratios=[1, 2], hspace=0.5)
  569. # Top: KDE
  570. ax_kde = fig.add_subplot(gs[0])
  572. # KDE plot data preparation
  573. body_results_dir = '/UK_BB/brainbody/bootstrap/body/results'
  574. file_impedance = 'body_composition_by_impedance_bootstrapped.csv'
  575. modality_impedance = 'body_composition_by_impedance'
  577. body_stacked_results_dir = '/UK_BB/brainbody/bootstrap/stack/scaled/results'
  578. file_stack = 'body_bootstrapped_scaled.csv'
  579. modality_stack = 'body'
  580. df_impedance = pd.read_csv(os.path.join(body_results_dir, file_impedance))
  581. df_stack = pd.read_csv(os.path.join(body_stacked_results_dir, file_stack))
  583. # Extract r values
  584. col_impedance = f'{modality_impedance}_r'
  585. col_stack = f'{modality_stack}_r'
  587. body_data = []
  588. if col_impedance in df_impedance.columns:
  589. r_vals = pd.to_numeric(df_impedance[col_impedance], errors='coerce').dropna()
  590. body_data.extend([{'modality': modality_impedance, 'r': r} for r in r_vals])
  592. if col_stack in df_stack.columns:
  593. r_vals_stack = pd.to_numeric(df_stack[col_stack], errors='coerce').dropna()
  594. body_data.extend([{'modality': modality_stack, 'r': r} for r in r_vals_stack])
  596. # Create DataFrame
  597. body_df = pd.DataFrame(body_data)
  599. # Compute medians
  600. median_impedance = body_df[body_df['modality'] == modality_impedance]['r'].median()
  601. median_stack = body_df[body_df['modality'] == modality_stack]['r'].median()
  603. # Plot KDE
  604. kde_palette = {
  605. 'body_composition_by_impedance': 'lightcoral',
  606. 'body': '#803342FF'
  607. }
  609. sns.kdeplot(
  610. data=body_df, x="r", hue="modality",
  611. fill=True, common_norm=False, palette=kde_palette,
  612. alpha=.5, linewidth=1, ax=ax_kde
  613. )
  615. # Get KDE heights at median positions
  616. r_impedance_vals = body_df[body_df['modality'] == modality_impedance]['r'].dropna()
  617. r_stack_vals = body_df[body_df['modality'] == modality_stack]['r'].dropna()
  619. kde_impedance = gaussian_kde(r_impedance_vals)
  620. kde_stack = gaussian_kde(r_stack_vals)
  622. # Evaluate KDE height at the median
  623. y_impedance = kde_impedance(median_impedance)[0]
  624. y_stack = kde_stack(median_stack)[0]
  626. # Draw short vertical lines within KDE height
  627. ax_kde.plot([median_impedance, median_impedance], [0, y_impedance],
  628. color='lightcoral', linestyle='--', linewidth=1)
  629. ax_kde.plot([median_stack, median_stack], [0, y_stack],
  630. color='#803342FF', linestyle='--', linewidth=1)
  632. # Get current y-axis upper limit
  633. y_max = ax_kde.get_ylim()[1]
  634. x_offset = 0.0
  635. ax_kde.text(median_impedance + x_offset, y_max * 0.99, f'{median_impedance:.2f}',
  636. color='black', fontsize=25, ha='center')
  637. ax_kde.text(median_stack + x_offset, y_max * 0.99, f'{median_stack:.2f}',
  638. color='black', fontsize=25, ha='center')
  640. sns.despine(left=True, ax=ax_kde)
  641. ax_kde.set_title("Bootstrap Distribution of Predictive Performance:\nTop Body Phenotype vs Composite Body Marker",
  642. fontsize=30, y=1.1, pad=20)
  644. # Add letter 'a' to top-left of KDE plot
  645. ax_kde.text(-0.1, 1.47, 'a', transform=ax_kde.transAxes,
  646. fontsize=40, va='top', ha='right') #fontweight='bold',
  649. ax_kde.set_xlabel("Pearson $r$", fontsize=27, labelpad=30)
  650. ax_kde.set_ylabel("Density", fontsize=25)
  652. # Adjust ticks
  653. r_min, r_max = body_df['r'].min(), body_df['r'].max()
  654. ax_kde.set_xticks(np.arange(round(r_min, 2), round(r_max + 0.01, 2), 0.01))
  655. ax_kde.tick_params(axis='x', labelsize=25)
  656. ax_kde.tick_params(axis='y', labelsize=25)
  658. # Remove legend
  659. ax_kde.legend_.remove()
  661. # Add labels under KDE distributions
  662. # Position labels at the bottom of the plot
  663. y_pos = -ax_kde.get_ylim()[1] * 0.15 # Position below x-axis
  665. # Add label for left distribution (body composition by impedance)
  666. ax_kde.text(median_impedance, y_pos, "Body composition by impedance", color ='lightcoral',
  667. fontsize=25, ha='center', va='top',
  668. transform=ax_kde.transData)
  670. # Add label for right distribution (composite body marker)
  671. ax_kde.text(median_stack, y_pos, "Composite body marker", color ='#803342FF',
  672. fontsize=25, ha='center', va='top',
  673. transform=ax_kde.transData)
  675. # Bottom: Violins
  676. # Create two subplots for the bottom row
  677. gs_bottom = gs[1].subgridspec(1, 2, width_ratios=[1, 1], wspace=0.15)
  678. ax1 = fig.add_subplot(gs_bottom[0]) # Left violin plot
  679. ax2 = fig.add_subplot(gs_bottom[1]) # Right violin plot
  681. # Styling patterns for violin plots
  682. violin_linewidth = 1
  683. violin_cut = 0
  684. violin_width = 0.8
  685. median_marker_size = 8
  686. median_color_left = '#A73933FF'
  687. median_color_right = '#53354DFF'
  688. violin_color_left = 'lightcoral'
  689. violin_color_right = '#DBC3D6FF'
  691. # Left: Individual body performance
  692. results_dir = '/UK_BB/brainbody/bootstrap/body/results'
  694. # Load files for left plot
  695. files = sorted([f for f in os.listdir(results_dir) if f.endswith('_bootstrapped.csv')])
  696. data_left = []
  698. for file in files:
  699. modality = file.replace('_bootstrapped.csv', '')
  700. df = pd.read_csv(os.path.join(results_dir, file))
  702. col = f'{modality}_r'
  703. if col in df.columns:
  704. r_vals = pd.to_numeric(df[col], errors='coerce').dropna()
  705. data_left.extend([{'modality': modality, 'r': r} for r in r_vals])
  706. else:
  707. print(f'Warning: expected column "{col}" not found in "{file}". Skipping.')
  709. df_left = pd.DataFrame(data_left)
  711. if not df_left.empty:
  712. # Compute median r per modality and sort descending
  713. medians_left = df_left.groupby('modality')['r'].median().sort_values(ascending=False)
  714. order_left = medians_left.index.tolist()
  716. # Determine display names for y-ticks via modality_map
  717. if 'modality_map' in globals() and isinstance(modality_map, dict):
  718. yticklabels_left = [modality_map.get(m, m) for m in order_left]
  719. else:
  720. yticklabels_left = order_left
  722. # Create left plot with consistent styling
  723. sns.violinplot(
  724. data=df_left,
  725. x='r', y='modality',
  726. order=order_left,
  727. ax=ax1,
  728. color=violin_color_left,
  729. linewidth=violin_linewidth,
  730. cut=violin_cut,
  731. inner=None,
  732. width=violin_width
  733. )
  735. # Add median markers for all modalities
  736. for j, modality in enumerate(order_left):
  737. median_val = medians_left[modality]
  738. ax1.scatter(median_val, j, color=median_color_left, s=median_marker_size, zorder=5)
  740. # Style left axes
  742. x_tick_fontsize = 25
  743. xlabel_fontsize = 27
  744. suptitle_fontsize = 28
  746. for spine in ax1.spines.values():
  747. spine.set_visible(True)
  749. # x ticks
  750. #ax1.xaxis.set_major_locator(MultipleLocator(0.1))
  751. ax1.set_xlim(left=0)
  752. ax1.tick_params(axis='both', which='both', length=4)
  753. ax1.tick_params(axis='x', labelsize=x_tick_fontsize)
  754. ax1.set_xlabel('Pearson $r$', fontsize=xlabel_fontsize, labelpad=10)
  756. x_min1, x_max1 = ax1.get_xlim()
  757. x_min1 = round(x_min1, 2)
  758. x_max1 = round(x_max1, 2)
  760. # Create fixed tick positions: min, min+0.1, min+0.2, max
  761. tick_positions = [0, x_min1 + 0.1, x_min1 + 0.2, x_min1 + 0.3, x_max1]
  762. # Set fixed locator
  763. from matplotlib.ticker import FixedLocator
  764. ax1.xaxis.set_major_locator(FixedLocator(tick_positions))
  766. # y ticks
  767. ax1.set_ylabel('')
  768. ax1.set_yticks(range(len(order_left)))
  769. ax1.set_yticklabels(yticklabels_left, fontsize=22)
  770. ax1.set_title('Individual\nBody Phenotype Performance', fontsize=suptitle_fontsize, pad=10)
  772. # Add letter 'b' to top-left of left violin plot
  773. ax1.text(-0.2, 1.15, 'b', transform=ax1.transAxes,
  774. fontsize=40, va='top', ha='right') #fontweight='bold',
  776. # Add CIs for left plot
  777. performance_bootstrapping_path = '/UK_BB/brainbody/bootstrap/body'
  778. performance_results_dir = os.path.join(performance_bootstrapping_path, 'results')
  779. ci_file = os.path.join(performance_results_dir, 'combined_bootstrap_CI_results.csv')
  780. if os.path.exists(ci_file):
  781. ci_df = pd.read_csv(ci_file)
  782. ci_df_r = ci_df[ci_df['Metric'] == 'r'] # Only get r values
  783. for j, modality in enumerate(order_left):
  784. if modality in ci_df_r['Modality'].values:
  785. ci_row = ci_df_r[ci_df_r['Modality'] == modality].iloc[0]
  786. lower_ci = ci_row['CI_lower']
  787. upper_ci = ci_row['CI_upper']
  789. # Add horizontal error bars at the bottom of each violin
  790. ax1.hlines(y=j + 0.2, xmin=lower_ci, xmax=upper_ci,
  791. color='black', linewidth=1, zorder=6)
  792. # Add vertical dashes for CI endpoints
  793. ax1.vlines([lower_ci, upper_ci], ymin=j + 0.1, ymax=j + 0.3,
  794. color='black', linewidth=1, zorder=6)
  796. # Clip violins to show only the TOP half for left plot
  797. xlim_left = ax1.get_xlim()
  798. violin_collections_left = [c for c in ax1.collections if isinstance(c, PolyCollection)]
  799. violin_collections_left = violin_collections_left[:len(order_left)]
  801. for j, coll in enumerate(violin_collections_left):
  802. rect = patches.Rectangle(
  803. (xlim_left[0], j - 0.4),
  804. xlim_left[1] - xlim_left[0],
  805. 0.4 + 1e-9,
  806. transform=ax1.transData,
  807. facecolor='none',
  808. edgecolor='none',
  809. linewidth=0
  810. )
  811. rect.set_visible(False)
  812. ax1.add_patch(rect)
  813. coll.set_clip_path(rect)
  815. # Right: Delta (Body vs Top Body)
  816. bootstrapping_path = '/UK_BB/brainbody/bootstrap/body'
  817. delta_results_dir = os.path.join(bootstrapping_path, 'delta_results')
  819. # Load delta results for right plot
  820. delta_dfs = {}
  821. if os.path.exists(delta_results_dir):
  822. delta_files = [f for f in os.listdir(delta_results_dir) if f.startswith('delta_') and f.endswith('.csv')]
  823. for f in delta_files:
  824. delta_dfs[f.replace('.csv', '')] = pd.read_csv(os.path.join(delta_results_dir, f))
  826. # Prepare data for right plot
  827. metric = 'delta_r'
  828. all_data_right = []
  830. for mod_name, df in delta_dfs.items():
  831. if metric in df.columns:
  832. comparison_name = mod_name.replace('delta_', '').replace('_vs_body_composition_by_impedance', '')
  833. comparison_name = modality_map.get(comparison_name, comparison_name)
  834. for val in df[metric]:
  835. all_data_right.append({'Comparison': comparison_name, 'Value': val})
  837. if all_data_right:
  838. plot_df = pd.DataFrame(all_data_right)
  839. sorted_order_right = plot_df.groupby('Comparison')['Value'].median().sort_values().index
  841. medians_right = plot_df.groupby('Comparison')['Value'].median().sort_values(ascending=False)
  843. # Create right plot
  844. sns.violinplot(
  845. data=plot_df,
  846. x='Value',
  847. y='Comparison',
  848. order=sorted_order_right,
  849. ax=ax2,
  850. color=violin_color_right,
  851. linewidth=violin_linewidth,
  852. cut=violin_cut,
  853. inner=None,
  854. width=violin_width
  855. )
  857. # Add median markers
  858. for j, modality in enumerate(sorted_order_right):
  859. median_val = medians_right[modality]
  860. ax2.scatter(median_val, j, color=median_color_right, s=median_marker_size, zorder=5)
  862. ax2.axvline(x=0, color='#A50026FF', linestyle='--', linewidth=0.7)
  863. ax2.set_ylabel('')
  865. # x ticks
  866. ax2.set_xlabel('Δ Pearson $r$', fontsize=xlabel_fontsize, labelpad=10)
  867. #ax2.xaxis.set_major_locator(MultipleLocator(0.1))
  868. #ax2.set_xlim(left=0)
  869. ax2.tick_params(axis='both', which='both', length=4)
  870. ax2.tick_params(axis='x', labelsize=x_tick_fontsize)
  872. x_min2, x_max2 = ax2.get_xlim()
  873. x_min2 = round(x_min2, 2)
  874. x_max2 = round(x_max2, 2)
  876. # Create fixed tick positions: min, min+0.1, min+0.2, max
  877. tick_positions = tick_positions = [x_min2, # show x_min
  878. x_min2 + 0.1,
  879. x_min2 + 0.2,
  880. x_min2 + 0.3,
  881. round(x_max2, 2)]
  883. # Set fixed locator
  884. from matplotlib.ticker import FixedLocator
  885. #ax2.xaxis.set_major_locator(FixedLocator(tick_positions))
  886. def format_tick(x, pos):
  887. if abs(x) < 1e-6:
  888. return "0.00"
  889. return f"{x:.2f}"
  891. from matplotlib.ticker import FixedLocator, FuncFormatter
  892. ax2.xaxis.set_major_formatter(FuncFormatter(format_tick))
  894. ax2.tick_params(axis='x', labelsize=x_tick_fontsize)
  895. ax2.set_title('Comparative\nPerformance', fontsize=suptitle_fontsize, pad=10)
  897. # Add CIs for right plot
  898. delta_ci_file = os.path.join(delta_results_dir, 'delta_CI_results.csv')
  899. if os.path.exists(delta_ci_file):
  900. delta_ci_df = pd.read_csv(delta_ci_file)
  901. delta_ci_df_r = delta_ci_df[delta_ci_df['Metric'] == 'r'].reset_index(drop=True)
  903. # Create a reverse mapping from display names to original names
  904. reverse_modality_map = {v: k for k, v in modality_map.items()}
  906. for j, modality in enumerate(sorted_order_right):
  907. original_modality_name = reverse_modality_map.get(modality, modality)
  908. ci_modality_name = f"{original_modality_name} vs body_composition_by_impedance"
  910. if ci_modality_name in delta_ci_df_r['Modality'].values:
  911. ci_row = delta_ci_df_r[delta_ci_df_r['Modality'] == ci_modality_name].iloc[0]
  912. lower_ci = ci_row['CI_lower']
  913. upper_ci = ci_row['CI_upper']
  915. # Add horizontal error bars at the bottom of each violin
  916. ax2.hlines(y=j - 0.2, xmin=lower_ci, xmax=upper_ci,
  917. color='black', linewidth=1, zorder=6)
  918. # Add vertical dashes for CI endpoints
  919. ax2.vlines([lower_ci, upper_ci], ymin=j - 0.3, ymax=j - 0.1,
  920. color='black', linewidth=1, zorder=6)
  921. else:
  922. print(f"NOT found: {ci_modality_name}")
  924. # Remove y-axis labels on the right plot but keep ticks
  925. ax2.set_yticks(range(len(sorted_order_right)))
  926. ax2.set_yticklabels([''] * len(sorted_order_right))
  928. # Control plot margins
  929. n_categories = len(sorted_order_right)
  930. ax2.set_ylim(-0.9, n_categories)
  932. # Clip violins for right plot
  933. xlim_right = ax2.get_xlim()
  934. violin_collections_right = [c for c in ax2.collections if isinstance(c, PolyCollection)]
  935. violin_collections_right = violin_collections_right[:len(sorted_order_right)]
  937. for j, coll in enumerate(violin_collections_right):
  938. rect = patches.Rectangle(
  939. (xlim_right[0], j),
  940. xlim_right[1] - xlim_right[0],
  941. 0.4 + 1e-9,
  942. transform=ax2.transData,
  943. facecolor='none',
  944. edgecolor='none',
  945. linewidth=0
  946. )
  947. rect.set_visible(False)
  948. ax2.add_patch(rect)
  949. coll.set_clip_path(rect)
  951. # Adjust y-axis limits for violin plots to be consistent
  952. max_categories = max(len(order_left), len(sorted_order_right))
  953. ax1.set_ylim(-0.5, max_categories - 0.5)
  954. ax2.set_ylim(-0.5, max_categories - 0.5)
  956. # Adjust y-ticks
  957. ax1.set_yticks(range(max_categories))
  958. ax2.set_yticks(range(max_categories))
  960. ax1.invert_yaxis()
  961. save_fig(['png', 'svg'], os.path.join(fig_path, 'final/Fig2'))
  962. plt.show()
  964. # %% [markdown]
  967. # %% [markdown]
  968. # # Figure 3
  970. # %% [markdown]
  971. # ## Figure 3a
  973. # %% [markdown]
  974. # Pearson *r* colorbar
  976. # %%
  977. # Normalize correlation values (0 to 1)
  978. norm = colors.Normalize(vmin=0, vmax=1)
  980. # Create a colormap that uses only the positive (blue) part of RdBu
  981. # RdBu goes from red (0.0) to white (0.5) to blue (1.0)
  982. positive_cmap = plt.colormaps['RdBu_r'].resampled(256)# Get the full colormap
  983. # Extract only the second half
  984. red_part = positive_cmap(np.linspace(0.5, 1, 256))
  985. positive_cmap = colors.LinearSegmentedColormap.from_list('positive_rdbu', red_part)
  987. # 1. Generate colorbar
  988. fig, ax = plt.subplots(figsize=(6, 1))
  989. fig.subplots_adjust(bottom=0.6)
  991. from matplotlib.colorbar import ColorbarBase
  992. cb1 = ColorbarBase(ax, cmap=positive_cmap, orientation='horizontal') #norm=norm,
  993. ticks = np.arange(0, 1.1, 0.2)
  994. cb1.set_ticks(ticks)
  995. cb1.set_ticklabels([f'{tick:.1f}' for tick in ticks])
  996. cb1.ax.tick_params(labelsize=25)
  997. cb1.set_label('Pearson $r$', fontsize=30)
  998. plt.savefig(os.path.join(fig_path, 'colorbar_horizontal.svg'), dpi=300)
  1000. # %% [markdown]
  1001. # Correlation dots
  1003. # %%
  1004. # Generate dot images
  1005. base_path = '/UK_BB/brainbody'
  1006. corr = pd.read_csv(os.path.join(base_path, 'feature_imp', 'feature_imp_body', 'combined', 'g_pred_from_body_mod_g_pred_stack_correlations.csv'))
  1007. for index, row in corr.iterrows():
  1008. modality = row['Modality']
  1009. r_value = row['Pearson r']
  1011. # Use positive colormap
  1012. color = positive_cmap(norm(r_value))
  1014. fig, ax = plt.subplots(figsize=(1, 1))
  1015. ax.add_patch(plt.Circle((0.5, 0.5), 0.4, color=color))
  1016. ax.set_xlim(0, 1)
  1017. ax.set_ylim(0, 1)
  1018. ax.axis('off')
  1020. # Clean filename
  1021. filename = modality.lower().replace(' ', '_').replace('&', 'and').replace('(', '').replace(')', '').replace(',', '').replace('/', '_')
  1022. filepath = os.path.join(fig_path, 'preps', 'body_performance', 'colorbar', 'positive', f'{filename}.png')
  1023. plt.savefig(filepath, dpi=300, bbox_inches='tight', pad_inches=0, transparent=True)
  1024. plt.close()
  1026. # %% [markdown]
  1027. # ## Figure 3b
  1029. # %% [markdown]
  1030. # Scatterplots of predicted vs observed g
  1032. # %%
  1033. # Scatterplots of predicted vs observed g: Config
  1034. bootstrapping_path = '/UK_BB/brainbody/bootstrap/stack'
  1035. modalities_to_plot = ['body', 'allmri', 'dwi', 'rs', 'smri', 'brain-body']
  1037. folds = range(0,5)
  1038. base_path = '/UK_BB/brainbody'
  1039. #base_path = 'Z:/IBu/UK_BB/brainbody'
  1041. dticmap = sns.light_palette("#0d648f", as_cmap=True)
  1042. rscmap = sns.light_palette("#4B6F5A99", as_cmap=True)
  1043. t1cmap = sns.light_palette("#f8b976", as_cmap=True)
  1044. stackcmap = sns.light_palette("#A73933FF", as_cmap=True)
  1045. bodycmap = sns.light_palette("#871C0FFF", as_cmap=True)
  1046. allmricmap = sns.light_palette("#6A659999", as_cmap=True)
  1048. modality_cmaps = {
  1049. 'dwi': dticmap,
  1050. 'rs': rscmap,
  1051. 'smri': t1cmap,
  1052. 'brain-body': stackcmap,
  1053. 'body': bodycmap,
  1054. 'allmri': allmricmap
  1055. }
  1057. modality_map = {
  1058. 'body': 'Composite Body',
  1059. 'allmri': 'Composite Brain',
  1060. 'dwi': 'dwMRI',
  1061. 'rs': 'rsMRI',
  1062. 'smri': 'sMRI',
  1063. 'brain-body': 'Whole-Body'
  1064. }
  1066. # -----------------------
  1067. # Load precomputed performance metrics AND bootstrap confidence intervals
  1068. brain_metrics = pd.read_excel(os.path.join(base_path, r'result/2level/2level_result-mean_brain_stack.xlsx'))
  1069. body_metrics = pd.read_excel(os.path.join(base_path, r'result/2level/2level_result-mean_body_stack_outer.xlsx'))
  1070. brain_body_metrics = pd.read_excel(os.path.join(base_path, r'result/2level/2level_result-mean_brain_and_body_stack_outer.xlsx'))
  1072. # Load bootstrap confidence intervals
  1073. bootstrap_ci_path = '/UK_BB/brainbody/bootstrap/stack/scaled/results/combined_bootstrap_CI_results_stacked_sorted_scaled.xlsx'
  1074. bootstrap_df = pd.read_excel(bootstrap_ci_path)
  1076. # Combine into one master dataframe
  1077. all_metrics_df = pd.concat([brain_metrics, body_metrics, brain_body_metrics], ignore_index=True)
  1079. # Create comprehensive annotation dictionary with bootstrap CIs
  1080. annotation_metrics = {}
  1082. for _, row in all_metrics_df.iterrows():
  1083. # Standardize modality names
  1084. modality_key = row['Modality'].lower()
  1086. # Replace specific patterns
  1087. modality_key = modality_key.replace('3 brain mri modalities', 'allmri')
  1088. modality_key = modality_key.replace('body physiology', 'body')
  1089. modality_key = modality_key.replace('brain mri and body physiology', 'brain-body')
  1090. modality_key = modality_key.replace('brain mri and body', 'brain-body')
  1091. modality_key = modality_key.replace('dwmri', 'dwi')
  1092. modality_key = modality_key.replace('rsmri', 'rs')
  1093. modality_key = modality_key.replace('smri', 'smri')
  1094. modality_key = modality_key.replace(' ', '-')
  1096. # Store the main metrics
  1097. annotation_metrics[modality_key] = {
  1098. 'r_value': row['Test Pearson r'],
  1099. 'r2_value': row['Test R2'],
  1100. 'r_ci_lower': None,
  1101. 'r_ci_upper': None,
  1102. 'r2_ci_lower': None,
  1103. 'r2_ci_upper': None
  1104. }
  1106. # Add bootstrap confidence intervals
  1107. for _, row in bootstrap_df.iterrows():
  1108. if row['Metric'] in ['R2', 'r']:
  1109. # Standardize modality names to match
  1110. modality_key = row['Modality'].lower()
  1111. modality_key = modality_key.replace('3 brain mri modalities', 'allmri')
  1112. modality_key = modality_key.replace('body physiology', 'body')
  1113. modality_key = modality_key.replace('brain mri and body physiology', 'brain-body')
  1114. modality_key = modality_key.replace('brain mri and body', 'brain-body')
  1115. modality_key = modality_key.replace('dwmri', 'dwi')
  1116. modality_key = modality_key.replace('rsmri', 'rs')
  1117. modality_key = modality_key.replace('smri', 'smri')
  1118. modality_key = modality_key.replace(' ', '-')
  1120. if modality_key in annotation_metrics:
  1121. if row['Metric'] == 'R2':
  1122. annotation_metrics[modality_key]['r2_ci_lower'] = row['CI_lower']
  1123. annotation_metrics[modality_key]['r2_ci_upper'] = row['CI_upper']
  1124. elif row['Metric'] == 'r':
  1125. annotation_metrics[modality_key]['r_ci_lower'] = row['CI_lower']
  1126. annotation_metrics[modality_key]['r_ci_upper'] = row['CI_upper']
  1128. print("Loaded annotation metrics for modalities:")
  1129. for modality, metrics in annotation_metrics.items():
  1130. print(f" {modality}: r={metrics['r_value']:.3f}, R²={metrics['r2_value']:.3f}")
  1131. if metrics['r_ci_lower'] is not None:
  1132. print(f" r CI: [{metrics['r_ci_lower']:.3f}, {metrics['r_ci_upper']:.3f}]")
  1133. if metrics['r2_ci_lower'] is not None:
  1134. print(f" R² CI: [{metrics['r2_ci_lower']:.3f}, {metrics['r2_ci_upper']:.3f}]")
  1136. # %%
  1137. # Plotting function
  1138. def plot_modality_scatter(
  1139. df, modality, save_dir=None,
  1140. cmap_name='viridis', edge_color='white',
  1141. marker_size=17, marker_alpha=0.6, line_lw=1.2, marker_linewidth=0.3,
  1142. title_fontsize=30, label_fontsize=35, tick_fontsize=30,
  1143. annotation_fontsize=None,
  1144. add_grid=False,
  1145. size=(6, 6)
  1146. ):
  1148. # Column names
  1149. x_col = 'g_pred_stack_test'
  1150. y_col = 'g_obs_test'
  1151. for col in (x_col, y_col):
  1152. if col not in df.columns:
  1153. raise ValueError(f"Required column '{col}' not found in data for {modality}. "
  1154. f"Available columns: {df.columns.tolist()}")
  1156. # Compute distance for color mapping (around the mean point)
  1157. x = pd.to_numeric(df[x_col], errors='coerce').astype(float)
  1158. y = pd.to_numeric(df[y_col], errors='coerce').astype(float)
  1159. x_mean, y_mean = x.mean(), y.mean()
  1160. dist_i = np.sqrt((y - y_mean)**2 + (x - x_mean)**2)
  1162. # Prepare figure/axes
  1163. fig, ax = plt.subplots(figsize=size)
  1165. # Scatter with colormap (use matplotlib scatter for 'c' + 'cmap' support)
  1166. sc = ax.scatter(
  1167. x, y, c=dist_i, cmap=cmap_name,
  1168. s=marker_size, alpha=marker_alpha,
  1169. edgecolors=edge_color, linewidth=marker_linewidth
  1170. )
  1172. # Global regression line
  1173. mask = np.isfinite(x) & np.isfinite(y)
  1174. if mask.sum() >= 2:
  1175. z = np.polyfit(x[mask], y[mask], 1)
  1176. p = np.poly1d(z)
  1177. # Sort x for a clean line
  1178. xs = np.linspace(x[mask].min(), x[mask].max(), 200)
  1179. ax.plot(xs, p(xs), alpha=0.8, linewidth=line_lw, color='#BE4A47FF', linestyle='-')
  1181. # Initialize annotation variables with default values
  1182. r_value = None
  1183. r2_value = None
  1184. r_ci_lower = None
  1185. r_ci_upper = None
  1186. r2_ci_lower = None
  1187. r2_ci_upper = None
  1188. annotation_text = "Default"
  1190. # Get annotation metrics from the pre-loaded dictionary
  1191. if modality in annotation_metrics:
  1192. metrics = annotation_metrics[modality]
  1193. r_value = metrics['r_value']
  1194. r2_value = metrics['r2_value']
  1195. r_ci_lower = metrics['r_ci_lower']
  1196. r_ci_upper = metrics['r_ci_upper']
  1197. r2_ci_lower = metrics['r2_ci_lower']
  1198. r2_ci_upper = metrics['r2_ci_upper']
  1200. # Create annotation text with confidence intervals
  1201. if r_ci_lower is not None and r2_ci_lower is not None:
  1202. annotation_text = (f'r = {r_value:.2f} [{r_ci_lower:.2f}, {r_ci_upper:.2f}]\n'
  1203. f'R² = {r2_value:.2f} [{r2_ci_lower:.2f}, {r2_ci_upper:.2f}]')
  1204. else:
  1205. # Fallback if CIs not available
  1206. annotation_text = f'r = {r_value:.2f}\nR² = {r2_value:.2f}'
  1207. else:
  1208. print(f"Warning: No metrics found for modality '{modality}'")
  1210. # Create annotation lines with proper error handling
  1211. annotation_lines = []
  1212. if r_value is not None and r_ci_lower is not None:
  1213. #annotation_lines.append(rf'$r$ = {r_value:.2f}, 95% CI [{r_ci_lower:.2f}, {r_ci_upper:.2f}]')
  1214. annotation_lines.append(rf'$r$ = {r_value:.2f}')
  1215. annotation_lines.append(rf'95% CI [{r_ci_lower:.2f}, {r_ci_upper:.2f}]')
  1216. elif r_value is not None:
  1217. annotation_lines.append(rf'$r$ = {r_value:.2f}')
  1219. # Aesthetics
  1220. sns.despine(top=True, right=True, ax=ax)
  1221. ax.set_xlabel('Predicted $g$-factor ($z$)', fontsize=label_fontsize)
  1222. ax.set_ylabel('$g$-factor derived from ESEM ($z$)', fontsize=label_fontsize)
  1223. ax.tick_params(axis='x', labelsize=tick_fontsize)
  1224. ax.tick_params(axis='y', labelsize=tick_fontsize)
  1226. # Title
  1227. display_name = modality_map.get(modality, modality) if 'modality_map' in globals() else modality
  1228. ax.set_title(display_name, fontsize=title_fontsize, y=1.14)
  1230. # Side annotation (top-left inside axes) - Only add if we have metrics
  1231. if annotation_lines:
  1232. ax.text(
  1233. 0.05, 0.9, "\n".join(annotation_lines),
  1234. transform=ax.transAxes, fontsize=annotation_fontsize,
  1235. va='bottom', ha='left'
  1236. )
  1238. # Major tick spacing
  1239. ax.xaxis.set_major_locator(MultipleLocator(1))
  1240. ax.yaxis.set_major_locator(MultipleLocator(1.5))
  1242. ax.set_xlim(-1.6, 1.5)
  1243. ax.set_ylim(-3.5, 3)
  1245. if add_grid:
  1246. ax.grid(True, alpha=0.25)
  1248. # Legend (only for lines)
  1249. handles, labels = ax.get_legend_handles_labels()
  1250. if handles:
  1251. ax.legend(loc='lower right', fontsize=12)
  1253. # Save and close
  1254. plt.tight_layout()
  1256. if save_dir is not None:
  1257. out_path = os.path.join(save_dir, f'scatter_{modality}_vs_observed.png')
  1258. plt.savefig(out_path, dpi=300, bbox_inches='tight')
  1260. plt.show()
  1261. plt.close()
  1263. # %%
  1264. # Plot scatterplots
  1265. modality_map = {
  1266. 'body': 'Composite Body',
  1267. 'allmri': 'Composite Brain',
  1268. 'dwi': 'dwMRI',
  1269. 'rs': 'rsMRI',
  1270. 'smri': 'sMRI',
  1271. 'brain-body': 'Whole-Body'
  1272. }
  1274. # Set up the save directory for plots
  1275. scatter_save_dir = '/UK_BB/brainbody/figures'
  1276. # Set annotation fontsize
  1277. annotation_fontsize = 25
  1279. # Process and plot each modality
  1280. for modality in modalities_to_plot:
  1281. print(f"\nGenerating scatter plot for {modality}...")
  1283. try:
  1284. # Load the combined test data for this modality
  1285. input_path = os.path.join(bootstrapping_path, 'scaled', f'{modality}_folds_combined_test.csv')
  1286. df = pd.read_csv(input_path)
  1288. print(f"Data shape for {modality}: {df.shape}")
  1289. print(f"Columns: {df.columns.tolist()}")
  1291. # Check if required columns exist
  1292. if 'g_obs_test' not in df.columns or 'g_pred_stack_test' not in df.columns:
  1293. print(f"Warning: Required columns not found for {modality}. Skipping.")
  1294. continue
  1296. # Get the appropriate colormap
  1297. cmap = modality_cmaps.get(modality, 'viridis')
  1299. # Create the scatter plot
  1300. plot_modality_scatter(
  1301. df=df,
  1302. modality=modality,
  1303. save_dir=scatter_save_dir,
  1304. cmap_name=cmap,
  1305. edge_color='black',
  1306. marker_size=50,
  1307. marker_alpha=0.9,
  1308. marker_linewidth=0.3,
  1309. line_lw=1.5,
  1310. title_fontsize=45,
  1311. label_fontsize=35,
  1312. tick_fontsize=35,
  1313. annotation_fontsize=35,
  1314. add_grid=False,
  1315. size=(8, 8) # Slightly larger for better visibility
  1316. )
  1318. print(f"Successfully plotted {modality}")
  1320. except Exception as e:
  1321. print(f"Error plotting {modality}: {str(e)}")
  1322. import traceback
  1323. traceback.print_exc()
  1325. print(f"\nAll scatter plots saved to: {scatter_save_dir}")
  1327. # %% [markdown]
  1328. # ## Figure 3c
  1330. # %% [markdown]
  1331. # Body feature importance
  1333. # %%
  1334. # Upload combined dataframe
  1335. base_path = '/UK_BB/brainbody'
  1336. output_path = os.path.join(base_path, 'feature_imp', 'feature_imp_body', 'combined')
  1337. combined = pd.read_csv(os.path.join(base_path, 'feature_imp/feature_imp_body/combined', 'body_features_g_pred_stack_combined.csv'))
  1338. corr_results = pd.read_csv(os.path.join(base_path, 'feature_imp/feature_imp_body/combined', 'body_features_g_pred_stack_correlations_detailed_sorted.csv'))
  1340. # %%
  1341. # 1. Filter correlations
  1342. significant_correlations = corr_results[
  1343. (corr_results['Pearson r'].abs() >= 0.198) &
  1344. (corr_results['p-value bonferroni'] < 0.05)
  1345. ].copy()
  1347. print(f"Found {len(significant_correlations)} features with |r| >= 0.198")
  1349. #significant_correlations = corr_results.copy()
  1351. # Check if the renamed phenotypes exist in the Phenotype column
  1352. if 'Body mass index (BMI) (Musculoskeletal)' in significant_correlations['Phenotype'].values:
  1353. print("Phenotype exists in results:")
  1354. print(significant_correlations[significant_correlations['Phenotype'] == 'Body mass index (BMI) (Musculoskeletal)'])
  1355. else:
  1356. print("Phenotype NOT found in results. Available phenotypes with 'BMI':")
  1357. bmi_phenotypes = significant_correlations[significant_correlations['Phenotype'].str.contains('BMI', na=False)]
  1358. for _, row in bmi_phenotypes.iterrows():
  1359. print(f" '{row['Phenotype']}': r = {row['Pearson r']:.3f}")
  1361. # 2. Identify which modality each feature belongs to
  1362. modality_mapping = {}
  1364. # Rename
  1365. modality_suffix_map = {
  1366. 'body_composition_by_impedance': ' (Impedance)',
  1367. 'body_composition_dxa': ' (DXA)',
  1368. 'musculoskeletal': ' (Musculoskeletal)',
  1369. 'cardiopulmonary': ' (Cardiopulmonary)',
  1370. 'arterial_stiffness': ' (Arterial Stiffness)'
  1371. }
  1373. duplicate_features_list = ['Trunk fat mass', 'Body mass index (BMI)', 'Weight', 'Pulse rate']
  1375. # Check one fold for each modality to get the feature names
  1376. for modality in modalities_body:
  1377. try:
  1378. # Load one fold's data for this modality
  1379. fold = 0 # Use fold 0 as reference
  1380. if modality == 'hearing':
  1381. scaled_path = os.path.join(
  1382. base_path,
  1383. 'hearing-vision',
  1384. 'folds',
  1385. f'fold_{fold}',
  1386. 'scaling',
  1387. f'{modality}_test_scaled_fold_{fold}.csv'
  1388. )
  1389. else:
  1390. scaled_path = os.path.join(
  1391. base_path,
  1392. 'body',
  1393. 'folds',
  1394. f'fold_{fold}',
  1395. 'scaling',
  1396. f'{modality}_test_scaled_fold_{fold}.csv'
  1397. )
  1399. # Load the modality features
  1400. modality_df = pd.read_csv(scaled_path)
  1401. modality_features = set(modality_df.columns)
  1403. # Apply renaming to features
  1404. if modality in modality_suffix_map:
  1405. suffix = modality_suffix_map[modality]
  1406. renamed_features = set()
  1407. for feature in modality_features:
  1408. if feature in duplicate_features_list:
  1409. renamed_features.add(f"{feature}{suffix}") #rename duplicate features
  1410. else:
  1411. renamed_features.add(feature)
  1412. modality_features = renamed_features
  1414. # Store the features for this modality
  1415. modality_mapping[modality] = modality_features
  1416. print(f"Phenotype '{modality}' has {len(modality_features)} features")
  1418. except Exception as e:
  1419. print(f"Error loading modality {modality}: {e}")
  1420. modality_mapping[modality] = set()
  1422. # Create reverse mapping for debugging
  1423. feature_to_modalities = {}
  1424. for modality, features in modality_mapping.items():
  1425. for feature in features:
  1426. if feature not in feature_to_modalities:
  1427. feature_to_modalities[feature] = []
  1428. feature_to_modalities[feature].append(modality)
  1430. # Check for duplicates (should be NONE now)
  1431. duplicate_features = {f: mods for f, mods in feature_to_modalities.items() if len(mods) > 1}
  1432. if duplicate_features:
  1433. print(f"\nFound {len(duplicate_features)} features in multiple modalities:")
  1434. for feature, mods in list(duplicate_features.items())[:5]:
  1435. print(f" '{feature}': {mods}")
  1436. else:
  1437. print(f"\nNo duplicate features found - all duplicates have been renamed!")
  1439. # 3 Check the distribution of modalities among significant features
  1440. def find_feature_modality(feature_name):
  1441. """Find which modality a feature belongs to"""
  1442. for modality, features in modality_mapping.items():
  1443. if feature_name in features:
  1444. return modality
  1445. return 'unknown' # If not found in any modality
  1447. significant_correlations['domain'] = significant_correlations['Phenotype'].apply(find_feature_modality)
  1449. # Check how many features map to each domain
  1450. print("\nDomain distribution:")
  1451. domain_counts = significant_correlations['domain'].value_counts()
  1452. print(domain_counts)
  1454. # Check for unknown modalities
  1455. unknown_features = significant_correlations[significant_correlations['domain'] == 'unknown']
  1456. if len(unknown_features) > 0:
  1457. print(f"\nWarning: {len(unknown_features)} features could not be mapped to a modality:")
  1458. for feature in unknown_features['Phenotype']:
  1459. print(f" - {feature}")
  1461. # Debug: Check why this feature isn't found
  1462. found_in = []
  1463. for modality, features in modality_mapping.items():
  1464. if feature in features:
  1465. found_in.append(modality)
  1466. if found_in:
  1467. print(f" Actually found in: {found_in}")
  1468. else:
  1469. print(f" Not found in any modality set")
  1470. else:
  1471. print("\nAll features successfully mapped to modalities")
  1473. # Check how many unique domains we have
  1474. unique_domains = significant_correlations['domain'].nunique()
  1475. print(f"\nUnique domains found: {unique_domains}")
  1476. print(f"Total features: {len(significant_correlations)}")
  1478. # 4 Add modality information to significant correlations
  1479. significant_correlations['domain'] = significant_correlations['Phenotype'].apply(find_feature_modality)
  1481. # Check for unknown modalities
  1482. unknown_features = significant_correlations[significant_correlations['domain'] == 'unknown']
  1483. if len(unknown_features) > 0:
  1484. print(f"\nWarning: {len(unknown_features)} features could not be mapped to a modality:")
  1485. for feature in unknown_features['Phenotype']:
  1486. print(f" - {feature}")
  1487. else:
  1488. print("\nAll features successfully mapped to modalities")
  1489. print(f"\nUnique domains found: {significant_correlations['domain'].nunique()}")
  1490. print(f"Total features: {len(significant_correlations)}")
  1493. # 5. Create the final mapping dataframe with sex-stratified data using Bonferroni correction
  1494. domain_corr_df = significant_correlations.rename(columns={
  1495. 'Phenotype': 'feature',
  1496. 'Pearson r': 'correlation',
  1497. 'p-value': 'p_value'
  1498. }).copy()
  1500. # Add sex-stratified columns (using Bonferroni-corrected p-values)
  1501. domain_corr_df['Pearson r male'] = significant_correlations['Pearson r male']
  1502. domain_corr_df['Pearson r female'] = significant_correlations['Pearson r female']
  1503. domain_corr_df['p_value_male'] = significant_correlations['p-value male bonferroni'] # Use Bonferroni-corrected
  1504. domain_corr_df['p_value_female'] = significant_correlations['p-value female bonferroni'] # Use Bonferroni-corrected
  1506. # Add significance flags for each sex using BONFERRONI CORRECTION (p < 0.05 after correction)
  1507. domain_corr_df['significant_male'] = domain_corr_df['p_value_male'] < 0.05
  1508. domain_corr_df['significant_female'] = domain_corr_df['p_value_female'] < 0.05
  1510. # Keep the overall significance (using Bonferroni-corrected)
  1511. domain_corr_df['p_value_overall_bonferroni'] = significant_correlations['p-value bonferroni']
  1512. domain_corr_df['significant'] = domain_corr_df['p_value_overall_bonferroni'] < 0.05
  1514. print(f"\nFinal domain mapping with Bonferroni-corrected sex-stratified data:")
  1515. print(domain_corr_df[['feature', 'correlation', 'domain', 'significant_male', 'significant_female']].head())
  1517. # Show summary of significant features
  1518. print(f"\nSignificant features summary:")
  1519. print(f"Significant in males: {domain_corr_df['significant_male'].sum()}")
  1520. print(f"Significant in females: {domain_corr_df['significant_female'].sum()}")
  1521. print(f"Significant in either sex: {(domain_corr_df['significant_male'] | domain_corr_df['significant_female']).sum()}")
  1523. # Save the mapping
  1524. mapping_output_path = os.path.join(base_path, 'feature_imp', 'feature_imp_body', 'combined', 'body_features_domain_mapping_sex_stratified_bonferroni.csv')
  1525. domain_corr_df.to_csv(mapping_output_path, index=False)
  1526. print(f"\nDomain mapping saved to: {mapping_output_path}")
  1528. # Show summary by domain
  1529. print("\nSummary by domain (Bonferroni-corrected):")
  1530. domain_summary = domain_corr_df.groupby('domain').agg({
  1531. 'feature': 'count',
  1532. 'correlation': ['mean', 'min', 'max'],
  1533. 'significant_male': 'sum',
  1534. 'significant_female': 'sum'
  1535. }).round(3)
  1536. print(domain_summary)
  1538. # 6 Create domain order
  1539. domain_order = [
  1540. 'body_composition_by_impedance',
  1541. 'bone_dxa',
  1542. 'abdominal_composition_mri_18_vars',
  1543. 'body_composition_dxa',
  1544. 'cardiopulmonary',
  1545. 'musculoskeletal',
  1546. 'abdominal_organ_composition_mri_13_vars',
  1547. 'metabolic',
  1548. 'renalhepatic',
  1549. 'pwa',
  1550. 'hearing',
  1551. 'carotid_ultrasound',
  1552. 'heart_mri',
  1553. 'kidneys_mri',
  1554. 'immune',
  1555. 'arterial_stiffness',
  1556. 'ecg_rest',
  1557. 'liver_mri',
  1558. 'bone_densitometry',
  1559. ]
  1560. #sorted(domain_corr_df['domain'].unique())
  1562. # Create a mapping of domain to y-position groups
  1563. domain_groups = {}
  1564. for i, domain in enumerate(domain_order):
  1565. domain_groups[domain] = i
  1567. # Assign y-position based on domain
  1568. domain_corr_df['domain_order'] = domain_corr_df['domain'].map(domain_groups)
  1570. # %%
  1571. # Rename features
  1572. feature_rename_map = {
  1573. 'L1-L4 average height': 'L1-L4 average height',
  1574. 'Forced expiratory volume (FEV1)': 'FEV1',
  1575. 'Sitting height': 'Sitting height',
  1576. 'Standing height': 'Standing height',
  1577. 'Seated height': 'Seated height',
  1578. 'Forced vital capacity (FVC)': 'FVC',
  1579. 'Height': 'Height',
  1580. 'Anterior thigh fat-free muscle volume (left)': 'Anterior thigh fat-free muscle vol (left)',
  1581. 'Femur wards BMD (bone mineral density) T-score (right)': 'Femur wards BMD T-score (right)',
  1582. 'Anterior thigh fat-free muscle volume (right)': 'Anterior thigh fat-free muscle vol (right)',
  1583. 'Head BMC (bone mineral content)': 'Head BMC',
  1584. 'IGF-1': 'IGF-1',
  1585. 'Peak expiratory flow (PEF)': 'PEF',
  1586. 'Head bone area': 'Head bone area',
  1587. 'Femur wards BMD (bone mineral density) T-score (left)': 'Femur wards BMD T-score (left)',
  1588. 'Femur wards BMD (bone mineral density) (right)': 'Femur wards BMD (right)',
  1589. 'Femur wards BMD (bone mineral density) (left)': 'Femur wards BMD (left)',
  1590. 'Femur neck BMD (bone mineral density) T-score (left)': 'Femur neck BMD T-score (left)',
  1591. 'Total BMD (bone mineral density) T-score': 'Total BMD T-score',
  1592. 'Femur neck BMD (bone mineral density) T-score (right)': 'Femur neck BMD T-score (right)',
  1593. 'Total thigh fat-free muscle volume': 'Total thigh fat-free muscle vol',
  1594. 'Pancreas volume': 'Pancreas volume',
  1595. 'Handgrip strength': 'Handgrip strength',
  1596. 'Femur neck BMD (bone mineral density) (left)': 'Femur neck BMD (left)',
  1597. 'Left kidney volume': 'Left kidney volume',
  1598. 'Femur neck BMD (bone mineral density) (right)': 'Femur neck BMD (right)',
  1599. 'Pelvis BMD (bone mineral density)': 'Pelvis BMD',
  1600. 'Posterior thigh fat-free muscle volume (right)': 'Posterior thigh fat-free muscle vol (right)',
  1601. 'Posterior thigh fat-free muscle volume (left)': 'Posterior thigh fat-free muscle vol (left)',
  1602. 'Gynoid lean mass': 'Gynoid lean mass',
  1603. 'Gynoid fat free mass': 'Gynoid fat free mass',
  1604. 'Mean arterial pressure during PWA': 'Mean arterial pressure during PWA',
  1605. 'Total trunk fat volume': 'Total trunk fat volume',
  1606. 'Trunk fat percentage': 'Trunk fat percentage',
  1607. 'Minimum carotid IMT (intima-medial thickness) at 120 degrees ': 'Minimum carotid IMT at 120°',
  1608. 'Minimum carotid IMT (intima-medial thickness) at 150 degrees ': 'Minimum carotid IMT at 150°',
  1609. 'Speech-reception-threshold (SRT) estimate (right)': 'SRT estimate (right)',
  1610. 'Total tissue fat percentage': 'Total tissue fat percentage',
  1611. 'Waist circumference': 'Waist circumference',
  1612. 'Speech-reception-threshold (SRT) estimate (left)': 'SRT estimate (left)',
  1613. 'Alkaline phosphatase': 'Alkaline phosphatase',
  1614. 'Stroke volume during PWA': 'Stroke volume during PWA',
  1615. 'Maximum carotid IMT (intima-medial thickness) at 210 degrees ': 'Maximum carotid IMT at 210°',
  1616. 'Mean carotid IMT (intima-medial thickness) at 210 degrees ': 'Mean carotid IMT at 210°',
  1617. 'Mean carotid IMT (intima-medial thickness) at 240 degrees ': 'Mean carotid IMT at 240°',
  1618. 'Maximum carotid IMT (intima-medial thickness) at 150 degrees ': 'Maximum carotid IMT at 150°',
  1619. 'Maximum carotid IMT (intima-medial thickness) at 240 degrees ': 'Maximum carotid IMT at 240°',
  1620. 'End systolic pressure during PWA': 'End systolic pressure during PWA',
  1621. 'Mean carotid IMT (intima-medial thickness) at 150 degrees ': 'Mean carotid IMT at 150°',
  1622. 'Maximum carotid IMT (intima-medial thickness) at 120 degrees ': 'Maximum carotid IMT at 120°',
  1623. 'Glycated haemoglobin (HbA1c)': 'Glycated haemoglobin (HbA1c)',
  1624. 'Mean carotid IMT (intima-medial thickness) at 120 degrees ': 'Mean carotid IMT at 120°',
  1625. 'Central augmentation pressure during PWA': 'Central augmentation pressure during PWA',
  1626. 'Cardiac index during PWA': 'Cardiac index during PWA',
  1627. 'Cardiac output during PWA': 'Cardiac output during PWA',
  1628. 'Trunk tissue fat percentage': 'Trunk tissue fat percentage',
  1629. 'Android tissue fat percentage': 'Android tissue fat percentage',
  1630. 'Visceral adipose tissue volume (VAT)': 'VAT vol',
  1631. 'Peripheral pulse pressure during PWA': 'Peripheral pulse pressure during PWA',
  1632. 'Systolic brachial blood pressure during PWA': 'Systolic brachial BP during PWA',
  1633. 'Systolic brachial blood pressure': 'Systolic brachial BP',
  1634. 'Visceral fat volume': 'Visceral fat vol',
  1635. 'Total abdominal adipose tissue index': 'Total abdominal adipose tissue index',
  1636. 'VAT (visceral adipose tissue) mass': 'VAT mass',
  1637. 'VAT (visceral adipose tissue) volume': 'VAT vol',
  1638. 'Central systolic blood pressure during PWA': 'Central systolic blood pressure during PWA',
  1639. 'Central pulse pressure during PWA': 'Central pulse pressure during PWA',
  1640. 'Weight-to-muscle ratio': 'Weight-to-muscle ratio',
  1641. 'Pancreas PDFF (fat fraction)': 'Pancreas PDFF (fat fraction)',
  1642. 'Cystatin C': 'Cystatin C',
  1643. 'Systolic blood pressure': 'Systolic BP',
  1644. 'Abdominal fat ratio': 'Abdominal fat ratio',
  1645. 'Posterior thigh muscle fat infiltration (MFI) (left)': 'Posterior thigh MFI (left)',
  1646. 'Posterior thigh muscle fat infiltration (MFI) (right)': 'Posterior thigh MFI (right)',
  1647. 'Anterior thigh muscle fat infiltration (MFI) (right)': 'Anterior thigh MFI (right)',
  1648. 'Anterior thigh muscle fat infiltration (MFI) (left)': 'Anterior thigh MFI (left)',
  1649. 'Muscle fat infiltration': 'Muscle fat infiltration',
  1650. }
  1652. # Apply renaming to Phenotype column
  1653. domain_corr_df['feature'] = domain_corr_df['feature'].map(feature_rename_map).fillna(domain_corr_df['feature'])
  1654. for i, (orig, new) in enumerate(zip(domain_corr_df['feature'], domain_corr_df['feature'])):
  1655. if orig != new:
  1656. print(f"'{orig}' -> '{new}'")
  1658. # %%
  1659. # Vertical feature importance plot
  1660. fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(60, 40),
  1661. gridspec_kw={'width_ratios': [1, 1]})
  1663. # Sort by domain and then by correlation
  1664. significant_df = domain_corr_df.sort_values(['domain_order', 'correlation'], ascending=[True, False])
  1665. significant_df = significant_df.reset_index(drop=True)
  1667. # Determine split point
  1668. total_features = len(significant_df)
  1669. split_idx = int(total_features / 2)
  1671. # Find the nearest domain boundary for a clean split
  1672. for i in range(split_idx, total_features):
  1673. if i == 0 or i == total_features:
  1674. continue
  1675. prev_domain = significant_df.iloc[i-1]['domain']
  1676. curr_domain = significant_df.iloc[i]['domain']
  1677. if curr_domain != prev_domain:
  1678. split_idx = i
  1679. print(f"Splitting at index {split_idx}, between domains '{prev_domain}' and '{curr_domain}'")
  1680. break
  1682. # Split the dataframe
  1683. df_left = significant_df.iloc[:split_idx].copy()
  1684. df_right = significant_df.iloc[split_idx:].copy()
  1686. print(f"Left plot: {len(df_left)} features, Right plot: {len(df_right)} features")
  1688. # Function to create a vertical bar plot
  1689. def create_vertical_subplot(ax, data_df, subplot_label=None, ylabel_fontsize=None):
  1690. """Create a vertical bar plot on given axes"""
  1692. # Create y positions with spacing between domains
  1693. y_positions = []
  1694. current_domain = None
  1695. spacing = 2.0
  1696. position = 0
  1698. for i in range(len(data_df)):
  1699. domain = data_df.iloc[i]['domain']
  1700. if domain != current_domain and current_domain is not None:
  1701. position += spacing
  1702. y_positions.append(position)
  1703. position += 1
  1704. current_domain = domain
  1706. y_positions = np.array(y_positions)
  1707. bar_height = 0.8
  1709. # Define color palette
  1710. custom_palette = [
  1711. "#A07756FF", '#BFCDD9FF', '#7F8C72FF', '#ADD0B5FF',
  1712. '#4A7169FF', '#BEB59CFF', '#4B644BFF', "#9B3F3FFF",
  1713. '#647D96FF', "#9970ABFF", '#226060FF', '#639CA4FF',
  1714. "#26588AFF", '#D2AD7CFF', "#2A775CFF","#C28C5FFF",
  1715. '#F4A582FF', "#B8858EFF", '#846D86FF'
  1716. ]
  1718. palette_dict = dict(zip(domain_corr_df['domain'].unique(), custom_palette))
  1720. # Create vertical bar plot
  1721. for i, (idx, row) in enumerate(data_df.iterrows()):
  1722. color = palette_dict[row['domain']]
  1723. ax.barh(y_positions[i], row['correlation'], height=bar_height,
  1724. color=color, alpha=0.9, edgecolor='black', linewidth=2.5)
  1726. # Add zero reference line
  1727. ax.axvline(x=0, color='grey', linewidth=3, linestyle='-', zorder=0)
  1729. # Set x-axis limits
  1730. max_corr = max(data_df['correlation'].max(), abs(data_df['correlation'].min()))
  1731. ax.set_xlim(-max_corr*1, max_corr*1)
  1732. ax.xaxis.set_major_locator(ticker.MultipleLocator(0.2))
  1734. # Manual domain offsets
  1735. domain_offset_overrides = {
  1736. 'body_composition_by_impedance': (-max_corr * 0.02, -1.6),
  1737. 'bone_dxa': (-max_corr * 0.02, -1.3),
  1738. 'abdominal_composition_mri_18_vars': (-max_corr * 0.02, -1.5),
  1739. 'body_composition_dxa': (-max_corr * 0.02, -1.4),
  1740. 'cardiopulmonary': (-max_corr * 0.02, -1.4),
  1741. 'musculoskeletal': (-max_corr * 0.02, -1.4),
  1742. 'abdominal_organ_composition_mri_13_vars': (-max_corr * 0.02, -1.4),
  1743. 'metabolic': (-max_corr * 0.02, -1.4),
  1744. 'renalhepatic': (-max_corr * 0.02, -1.4),
  1745. 'pwa': (-max_corr * 0.02, -1.4),
  1746. 'hearing': (-max_corr * 0.02, -1.4),
  1747. 'carotid_ultrasound': (-max_corr * 0.02, -1.4)
  1748. }
  1750. # Initialize domain_offsets dictionary
  1751. domain_offsets = {}
  1752. current_domain = None
  1753. domain_start_idx = 0
  1755. # Process each row to calculate offsets
  1756. for i in range(len(data_df)):
  1757. domain = data_df.iloc[i]['domain']
  1758. if domain != current_domain:
  1759. if current_domain is not None:
  1760. # Use manual override if defined, otherwise use default
  1761. if current_domain in domain_offset_overrides:
  1762. x_offset, y_offset = domain_offset_overrides[current_domain]
  1763. else:
  1764. # Default offsets for domains not in overrides
  1765. x_offset = -max_corr * 0.15 # Default horizontal
  1766. y_offset = -0.2 # Default vertical
  1768. # Store the calculated offset
  1769. domain_offsets[current_domain] = (x_offset, y_offset)
  1771. domain_start_idx = i
  1772. current_domain = domain
  1774. # Handle the last domain
  1775. if current_domain is not None:
  1776. if current_domain in domain_offset_overrides:
  1777. x_offset, y_offset = domain_offset_overrides[current_domain]
  1778. else:
  1779. # Default offsets for domains not in overrides
  1780. x_offset = -max_corr * 0.15
  1781. y_offset = -0.2
  1783. domain_offsets[current_domain] = (x_offset, y_offset)
  1785. # Add domain annotations using the calculated offsets
  1786. current_domain = None
  1787. domain_start_idx = 0
  1789. for i in range(len(data_df)):
  1790. domain = data_df.iloc[i]['domain']
  1791. if domain != current_domain:
  1792. if current_domain is not None:
  1793. x_offset, y_offset = domain_offsets.get(current_domain, (-max_corr * 0.15, -0.2))
  1794. annotation_y = y_positions[domain_start_idx] + y_offset
  1796. ax.text(x_offset, annotation_y, modality_map[current_domain],
  1797. ha='right', va='center', fontsize=65, fontweight='bold', rotation=0,
  1798. bbox=dict(boxstyle="round,pad=0.1", facecolor=palette_dict[current_domain], alpha=0.5, linewidth=3))
  1800. domain_start_idx = i
  1801. current_domain = domain
  1803. # Add annotation for the very last domain
  1804. if current_domain is not None:
  1805. x_offset, y_offset = domain_offsets.get(current_domain, (-max_corr * 0.15, -0.2))
  1806. annotation_y = y_positions[domain_start_idx] + y_offset
  1808. ax.text(x_offset, annotation_y, modality_map[current_domain],
  1809. ha='right', va='center', fontsize=65, fontweight='bold', rotation=0,
  1810. bbox=dict(boxstyle="round,pad=0.1", facecolor=palette_dict[current_domain], alpha=0.5))
  1812. # Add separation lines between domains
  1813. current_domain = None
  1814. for i in range(len(data_df)):
  1815. domain = data_df.iloc[i]['domain']
  1816. if domain != current_domain:
  1817. if current_domain is not None:
  1818. gap_center = (y_positions[i-1] + y_positions[i]) / 2
  1819. ax.axhline(y=gap_center, color='black', xmin=0.55, xmax=1.0, linestyle=':', alpha=0.7, linewidth=2)
  1820. current_domain = domain
  1822. # Set y-axis limits and labels
  1823. ax.set_ylim(y_positions[0] - 1, y_positions[-1] + 1)
  1824. ax.set_yticks(y_positions)
  1825. ax.set_yticklabels(data_df['feature'].tolist(), fontsize=ylabel_fontsize)
  1827. # Customize axes
  1828. ax.tick_params(axis='x',
  1829. labelsize=80,
  1830. bottom=True, # bottom ticks
  1831. length=10,
  1832. width=2,
  1833. color='black')
  1834. ax.set_xlabel("Pearson $r$", fontsize=90, labelpad=20)
  1836. # Customize spines
  1837. ax.spines['left'].set_linewidth(0)
  1838. ax.spines['bottom'].set_linewidth(2)
  1839. ax.spines['top'].set_linewidth(0)
  1840. ax.spines['right'].set_linewidth(0)
  1842. # Invert y-axis
  1843. ax.invert_yaxis()
  1845. # Add subplot label (A, B)
  1846. if subplot_label:
  1847. ax.text(0.02, 0.98, subplot_label, transform=ax.transAxes,
  1848. fontsize=44, fontweight='bold', va='top',
  1849. bbox=dict(boxstyle="round,pad=0.3", facecolor='white', alpha=0.7))
  1851. return max_corr
  1853. # Create left subplot
  1854. print("Creating left subplot...")
  1855. max_corr_left = create_vertical_subplot(ax1, df_left, subplot_label=None, ylabel_fontsize=55)
  1857. # Create right subplot
  1858. print("Creating right subplot...")
  1859. max_corr_right = create_vertical_subplot(ax2, df_right, subplot_label=None, ylabel_fontsize=55)
  1861. # Ensure both plots have consistent x-axis limits for comparison
  1862. max_corr_both = max(max_corr_left, max_corr_right)
  1863. ax1.set_xlim(-max_corr_both*1.3, max_corr_both*1.1)
  1864. ax2.set_xlim(-max_corr_both*1.3, max_corr_both*1.1)
  1867. # Share y-axis scale
  1868. # ax1.sharey(ax2)
  1870. fig.suptitle('Feature Importance:\nAssociations between Composite Body Marker\nand Each Body Physiology Measure (|$r$| ≥ 0.2)',
  1871. fontsize=150,
  1872. #fontweight='bold',
  1873. y=1.05,
  1874. x=0.57,
  1875. verticalalignment='top')
  1877. # Adjust tight_layout to make space for title
  1878. plt.tight_layout(rect=[0, 0, 1, 0.96]) # Leave top 4% for title
  1880. # Adjust layout
  1881. plt.tight_layout()
  1883. # Save the figure
  1884. plt.savefig(os.path.join(fig_path, 'Fig3c_body_features_g_pred_stack_annotated_split_vertical.svg'),
  1885. bbox_inches="tight",
  1886. pad_inches=1,
  1887. transparent=False,
  1888. facecolor="w",
  1889. edgecolor='w',
  1890. format='svg',
  1891. dpi=300)
  1893. plt.show()
  1894. print(f"Figure saved with {len(df_left)} features in left plot and {len(df_right)} features in right plot")
  1896. # %% [markdown]
  1897. # # Figure 4
  1899. # %% [markdown]
  1900. # ## Figure 4a
  1902. # %% [markdown]
  1903. # Composite body performance versus brain and whole-body markers
  1905. # %%
  1906. # Prepare data for KDE plot of delta_r values
  1907. import warnings
  1908. warnings.simplefilter(action='ignore', category=FutureWarning)
  1910. bootstrapping_path = '/UK_BB/brainbody/bootstrap/body_vs_all'
  1911. delta_dfs = {}
  1912. delta_dir_stacked = os.path.join(bootstrapping_path, 'delta_results_stacked')
  1913. if os.path.exists(delta_dir_stacked):
  1914. delta_files = [f for f in os.listdir(delta_dir_stacked) if f.startswith('delta_') and f.endswith('.csv')]
  1915. for f in delta_files:
  1916. delta_dfs[f.replace('.csv', '')] = pd.read_csv(os.path.join(delta_dir_stacked, f))
  1919. kde_palette = {
  1920. 'dwMRI': '#0d648f',
  1921. 'rsMRI': '#4B6F5A99',
  1922. 'sMRI': '#f8b976',
  1923. '3 Brain MRI Modalities': '#6A659999',
  1924. '3 Brain MRI Modalities and Body': '#A73933FF'
  1925. }
  1927. combined_data = []
  1928. delta_modality_mapping = {
  1929. 'delta_dwi_vs_body': 'dwMRI',
  1930. 'delta_rs_vs_body': 'rsMRI',
  1931. 'delta_smri_vs_body': 'sMRI',
  1932. 'delta_allmri_vs_body': '3 Brain MRI Modalities',
  1933. 'delta_brain-body_vs_body': '3 Brain MRI Modalities and Body'
  1934. }
  1936. for delta_name, df in delta_dfs.items():
  1937. if 'delta_r' in df.columns and delta_name in delta_modality_mapping:
  1938. display_name = delta_modality_mapping[delta_name]
  1939. for delta_r_val in df['delta_r']:
  1940. combined_data.append({'modality': display_name, 'delta_r': delta_r_val})
  1942. combined_df = pd.DataFrame(combined_data)
  1944. # %%
  1945. # KDE plot: Composite body vs brain and whole-body markers
  1946. fig = plt.figure(figsize=(15,3))
  1948. kde = sns.kdeplot(
  1949. data=combined_df, x="delta_r", hue="modality",
  1950. fill=True, common_norm=False, palette=kde_palette,
  1951. alpha=.5, linewidth=1, legend=False
  1952. )
  1954. # Add median lines and labels
  1955. modalities = combined_df['modality'].unique()
  1956. for mod in modalities:
  1957. median_delta = combined_df[combined_df['modality'] == mod]['delta_r'].median()
  1958. #plt.axvline(median_delta, color='grey', linestyle='--', linewidth=0.7)
  1960. # Different x-offsets to prevent overlap
  1961. if median_delta < 0.01:
  1962. x_offset = -0.008
  1963. ha = 'left'
  1964. elif median_delta < 0.05:
  1965. x_offset = -0.009
  1966. ha = 'left'
  1967. elif median_delta < 0.065:
  1968. x_offset = -0.009
  1969. ha = 'left'
  1970. elif median_delta > 0.07:
  1971. x_offset = 0.004
  1972. ha = 'left'
  1973. else: # Larger values
  1974. x_offset = 0.005
  1975. ha = 'left'
  1977. plt.text(median_delta + x_offset, plt.gca().get_ylim()[1] * 0.88 + 8,
  1978. f'{median_delta:.3f}', color='black', fontsize=25, ha=ha)
  1980. # Add zero reference line
  1981. plt.axvline(x=0, color='#32363FFF', linestyle='--', linewidth=2, alpha=1) #BE4A47FF
  1982. # Add short vertical lines at KDE peak for each modality
  1983. for mod in modalities:
  1984. r_vals = combined_df[combined_df['modality'] == mod]['delta_r'].dropna()
  1985. median_delta = r_vals.median()
  1987. # Estimate KDE height at median
  1988. kde = gaussian_kde(r_vals)
  1989. y_peak = kde(median_delta)[0]
  1991. # Draw vertical line from y=0 to KDE peak
  1992. plt.plot([median_delta, median_delta], [0, y_peak], color=kde_palette.get(mod, 'grey'), linestyle='--', linewidth=1)
  1994. # ------------
  1995. sns.despine(left=True)
  1996. plt.title("Comparative Performance of Stacked Models", fontsize=30, y=1.3) #Body Stacked vs Brain Stacked
  1997. plt.xlabel("Δ Pearson $r$", fontsize=30, labelpad=50)
  1998. plt.ylabel("Density", fontsize=25)
  1999. plt.yticks(fontsize=20)
  2000. #plt.xticks(fontsize=20)
  2002. # Adjust x-ticks for delta values
  2003. x_min = combined_df['delta_r'].min()
  2004. x_max = combined_df['delta_r'].max()
  2005. x_ticks = np.arange(round(x_min, 2), round(x_max, 2) + 0.02, 0.02)
  2006. plt.xticks(ticks=x_ticks, labels=[f'{x:.2f}' for x in x_ticks], fontsize=20)
  2008. plt.savefig(os.path.join(fig_path, 'Fig4a_delta_r_brain_vs_body_stacked_scaled_KDE.svg'),
  2009. bbox_inches="tight",
  2010. #pad_inches=1,
  2011. transparent=False,
  2012. facecolor="w",
  2013. edgecolor='w',
  2014. orientation='landscape',
  2015. format='svg')
  2017. plt.show()
  2019. # %% [markdown]
  2020. # ## Figure 4b
  2022. # %% [markdown]
  2023. # Combined scatterplots of predicted vs observed g
  2025. # %%
  2026. # Combined scatterplots of predicted vs observed g: Config
  2027. bootstrapping_path = '/UK_BB/brainbody/bootstrap/stack'
  2028. modalities_to_plot = ['body', 'allmri', 'dwi', 'rs', 'smri', 'brain-body']
  2030. folds = range(0,5)
  2031. base_path = '/UK_BB/brainbody'
  2033. dticmap = sns.light_palette("#0d648f", as_cmap=True)
  2034. rscmap = sns.light_palette("#4B6F5A99", as_cmap=True)
  2035. t1cmap = sns.light_palette("#f8b976", as_cmap=True)
  2036. stackcmap = sns.light_palette("#A73933FF", as_cmap=True)
  2037. bodycmap = sns.light_palette("#871C0FFF", as_cmap=True)
  2038. allmricmap = sns.light_palette("#6A659999", as_cmap=True)
  2040. modality_cmaps = {
  2041. 'dwi': dticmap,
  2042. 'rs': rscmap,
  2043. 'smri': t1cmap,
  2044. 'brain-body': stackcmap,
  2045. 'body': bodycmap,
  2046. 'allmri': allmricmap
  2047. }
  2049. modality_map = {
  2050. 'body': 'Composite Body',
  2051. 'allmri': 'Composite Brain',
  2052. 'dwi': 'dwMRI',
  2053. 'rs': 'rsMRI',
  2054. 'smri': 'sMRI',
  2055. 'brain-body': 'Body+Brain'
  2056. }
  2059. # -----------------------
  2060. # Load precomputed performance metrics AND bootstrap confidence intervals
  2061. brain_metrics = pd.read_excel(os.path.join(base_path, r'result/2level/2level_result-mean_brain_stack.xlsx'))
  2062. body_metrics = pd.read_excel(os.path.join(base_path, r'result/2level/2level_result-mean_body_stack_outer.xlsx'))
  2063. brain_body_metrics = pd.read_excel(os.path.join(base_path, r'result/2level/2level_result-mean_brain_and_body_stack_outer.xlsx'))
  2065. # Load bootstrap confidence intervals
  2066. bootstrap_ci_path = '/UK_BB/brainbody/bootstrap/stack/scaled/results/combined_bootstrap_CI_results_stacked_sorted_scaled.xlsx'
  2067. bootstrap_df = pd.read_excel(bootstrap_ci_path)
  2069. # Combine into one master dataframe
  2070. all_metrics_df = pd.concat([brain_metrics, body_metrics, brain_body_metrics], ignore_index=True)
  2072. # Create comprehensive annotation dictionary with bootstrap CIs
  2073. annotation_metrics = {}
  2075. for _, row in all_metrics_df.iterrows():
  2076. # Standardize modality names
  2077. modality_key = row['Modality'].lower()
  2079. # Replace specific patterns
  2080. modality_key = modality_key.replace('3 brain mri modalities', 'allmri')
  2081. modality_key = modality_key.replace('body physiology', 'body')
  2082. modality_key = modality_key.replace('brain mri and body physiology', 'brain-body')
  2083. modality_key = modality_key.replace('brain mri and body', 'brain-body')
  2084. modality_key = modality_key.replace('dwmri', 'dwi')
  2085. modality_key = modality_key.replace('rsmri', 'rs')
  2086. modality_key = modality_key.replace('smri', 'smri')
  2087. modality_key = modality_key.replace(' ', '-')
  2089. # Store the main metrics
  2090. annotation_metrics[modality_key] = {
  2091. 'r_value': row['Test Pearson r'],
  2092. 'r2_value': row['Test R2'],
  2093. 'r_ci_lower': None,
  2094. 'r_ci_upper': None,
  2095. 'r2_ci_lower': None,
  2096. 'r2_ci_upper': None
  2097. }
  2099. # Add bootstrap confidence intervals
  2100. for _, row in bootstrap_df.iterrows():
  2101. if row['Metric'] in ['R2', 'r']:
  2102. # Standardize modality names to match
  2103. modality_key = row['Modality'].lower()
  2104. modality_key = modality_key.replace('3 brain mri modalities', 'allmri')
  2105. modality_key = modality_key.replace('body physiology', 'body')
  2106. modality_key = modality_key.replace('brain mri and body physiology', 'brain-body')
  2107. modality_key = modality_key.replace('brain mri and body', 'brain-body')
  2108. modality_key = modality_key.replace('dwmri', 'dwi')
  2109. modality_key = modality_key.replace('rsmri', 'rs')
  2110. modality_key = modality_key.replace('smri', 'smri')
  2111. modality_key = modality_key.replace(' ', '-')
  2113. if modality_key in annotation_metrics:
  2114. if row['Metric'] == 'R2':
  2115. annotation_metrics[modality_key]['r2_ci_lower'] = row['CI_lower']
  2116. annotation_metrics[modality_key]['r2_ci_upper'] = row['CI_upper']
  2117. elif row['Metric'] == 'r':
  2118. annotation_metrics[modality_key]['r_ci_lower'] = row['CI_lower']
  2119. annotation_metrics[modality_key]['r_ci_upper'] = row['CI_upper']
  2121. print("Loaded annotation metrics for modalities:")
  2122. for modality, metrics in annotation_metrics.items():
  2123. print(f" {modality}: r={metrics['r_value']:.3f}, R²={metrics['r2_value']:.3f}")
  2124. if metrics['r_ci_lower'] is not None:
  2125. print(f" r CI: [{metrics['r_ci_lower']:.3f}, {metrics['r_ci_upper']:.3f}]")
  2126. if metrics['r2_ci_lower'] is not None:
  2127. print(f" R² CI: [{metrics['r2_ci_lower']:.3f}, {metrics['r2_ci_upper']:.3f}]")
  2129. # %%
  2130. # Combined scatterplots
  2131. ordered_modalities = ['smri', 'rs', 'dwi', 'brain-body', 'allmri']
  2133. # Create figure with 1 row and 5 columns
  2134. fig, axes = plt.subplots(1, 5, figsize=(35, 8))
  2136. for idx, modality in enumerate(ordered_modalities):
  2137. if idx >= len(axes):
  2138. break
  2140. try:
  2141. # Load data
  2142. input_path = os.path.join(bootstrapping_path, f'{modality}_folds_combined_test.csv')
  2143. df = pd.read_csv(input_path)
  2145. # Get colormap
  2146. cmap = modality_cmaps[modality]
  2148. # Plot on subplot
  2149. ax = axes[idx]
  2150. x_col = 'g_pred_stack_test'
  2151. y_col = 'g_obs_test'
  2153. x = pd.to_numeric(df[x_col], errors='coerce').astype(float)
  2154. y = pd.to_numeric(df[y_col], errors='coerce').astype(float)
  2156. # Compute distance for coloring
  2157. x_mean, y_mean = x.mean(), y.mean()
  2158. dist_i = np.sqrt((y - y_mean)**2 + (x - x_mean)**2)
  2160. # Scatter plot - updated parameters
  2161. sc = ax.scatter(
  2162. x, y, c=dist_i, cmap=cmap,
  2163. s=50, alpha=0.9,
  2164. edgecolors='black', linewidth=0.3
  2165. )
  2167. # Regression line
  2168. mask = np.isfinite(x) & np.isfinite(y)
  2169. if mask.sum() >= 2:
  2170. z = np.polyfit(x[mask], y[mask], 1)
  2171. p = np.poly1d(z)
  2172. xs = np.linspace(x[mask].min(), x[mask].max(), 200)
  2173. ax.plot(xs, p(xs), alpha=0.8, linewidth=2, color='#BE4A47FF', linestyle='-')
  2175. # Set axis limits and ticks
  2176. ax.set_xlim(-3, 2)
  2177. ax.set_ylim(-3.5, 3)
  2178. ax.xaxis.set_major_locator(MultipleLocator(1))
  2179. ax.yaxis.set_major_locator(MultipleLocator(1))
  2181. # Add annotation with updated format
  2182. if modality in annotation_metrics:
  2183. metrics = annotation_metrics[modality]
  2184. r_value = metrics['r_value']
  2185. r_ci_lower = metrics.get('r_ci_lower')
  2186. r_ci_upper = metrics.get('r_ci_upper')
  2188. # Updated annotation text format
  2189. if r_ci_lower is not None:
  2190. annotation_text = f'$r$ = {r_value:.2f}'
  2191. annotation_text += f'\n95% CI [{r_ci_lower:.2f}, {r_ci_upper:.2f}]'
  2192. else:
  2193. annotation_text = rf'$r$ = {r_value:.2f}'
  2195. ax.text(
  2196. 0.05, 0.91, annotation_text,
  2197. transform=ax.transAxes, fontsize=38,
  2198. va='bottom', ha='left'
  2199. )
  2201. # Set labels with updated format
  2202. if idx == 0: # Only for first plot
  2203. ax.set_ylabel('$g$-factor\nderived from ESEM ($z$)', fontsize=50)
  2204. else:
  2205. ax.set_ylabel('') # Empty for plots 2-5
  2206. ax.set_yticklabels([])
  2208. # Set title with bbox edge color matching modality cmap
  2209. title_color = cmap(0.5)
  2210. ax.set_title(modality_map.get(modality, modality), fontsize=50, y=1.2,
  2211. bbox=dict(facecolor='white', alpha=0.9, edgecolor=title_color,
  2212. linewidth=2, boxstyle='round,pad=0.3'))
  2214. # Remove frame (spines) - updated formatting
  2215. sns.despine(top=True, right=True, ax=ax)
  2217. # Set tick label sizes - updated formatting
  2218. ax.tick_params(axis='x', labelsize=40)
  2219. ax.tick_params(axis='y', labelsize=40)
  2221. except Exception as e:
  2222. print(f"Error in combined plot for {modality}: {str(e)}")
  2223. axes[idx].text(0.5, 0.5, f"Error\n{modality}", ha='center', va='center', transform=axes[idx].transAxes)
  2225. # Remove empty subplots if any
  2226. for idx in range(len(ordered_modalities), len(axes)):
  2227. fig.delaxes(axes[idx])
  2229. # Add global x-label
  2230. fig.text(0.5, -0.12, 'Predicted $g$-factor ($z$)', ha='center', fontsize=50)
  2232. plt.tight_layout()
  2234. scatter_save_dir = '/UK_BB/brainbody/figures'
  2235. plt.savefig(os.path.join(scatter_save_dir, 'Fig4b_all_modalities_combined.png'), dpi=300, bbox_inches='tight', pad_inches=0.5)
  2236. plt.show()
  2238. # %% [markdown]
  2239. # ## Figure 4c
  2241. # %% [markdown]
  2242. # Feature importance: heatmaps dwMRI + rsMRI
  2244. # %%
  2245. # Functions to transform 1D correlations dataframe to 2D correlation matrix
  2246. def long_to_square_matrix_dwi(correlation, coordinates):
  2247. regions = coordinates['name'].values
  2249. n = len(regions)
  2250. corr_matrix = np.zeros((n, n))
  2251. pval_matrix = np.ones((n, n))
  2253. name_to_idx = {name: idx for idx, name in enumerate(regions)}
  2255. for _, row in correlation.iterrows():
  2256. i = name_to_idx[row['region_i']]
  2257. j = name_to_idx[row['region_j']]
  2258. corr_matrix[i, j] = row['Pearson r']
  2259. corr_matrix[j, i] = row['Pearson r']
  2260. pval_matrix[i, j] = row['p-value']
  2261. pval_matrix[j, i] = row['p-value']
  2263. return (
  2264. pd.DataFrame(corr_matrix, index=regions, columns=regions),
  2265. pd.DataFrame(pval_matrix, index=regions, columns=regions)
  2266. )
  2268. def long_to_square_matrix_rs(correlation, coordinates):
  2269. regions = coordinates['name'].values
  2270. n = len(regions)
  2272. corr_matrix = np.full((n, n), np.nan)
  2273. pval_matrix = np.full((n, n), np.nan)
  2275. np.fill_diagonal(corr_matrix, 1.0)
  2276. np.fill_diagonal(pval_matrix, 0.0)
  2278. name_to_idx = {name: idx for idx, name in enumerate(regions)}
  2280. for _, row in correlation.iterrows():
  2281. i = name_to_idx[row['region_i']]
  2282. j = name_to_idx[row['region_j']]
  2283. corr_matrix[i, j] = row['Pearson r']
  2284. corr_matrix[j, i] = row['Pearson r']
  2285. pval_matrix[i, j] = row['p-value']
  2286. pval_matrix[j, i] = row['p-value']
  2288. return (
  2289. pd.DataFrame(corr_matrix, index=regions, columns=regions),
  2290. pd.DataFrame(pval_matrix, index=regions, columns=regions)
  2291. )
  2293. # Load data
  2294. base_path = '/UK_BB/brainbody'
  2295. schaefer_corr_dwi = '/UK_BB/brainbody/feature_imp/feature_imp_brain/dwmri/Schaefer7n200p_Tian_S1_Streamline_Count_i2_corr_with_g_stack_with_regions.csv'
  2296. schaefer_corr_rs = '/UK_BB/brainbody/feature_imp/feature_imp_brain/rsmri/Schaefer7n200p_Tian_S1_Full_Corr_i2_corr_with_g_stack_with_regions.csv'
  2298. corr_dwi = pd.read_csv(schaefer_corr_dwi)
  2299. corr_rs = pd.read_csv(schaefer_corr_rs)
  2300. coords = pd.read_csv('/UK_BB/brainbody/feature_imp/feature_imp_brain/mni_coords.csv')
  2302. # Read region order
  2303. with open(os.path.join(base_path, 'feature_imp/feature_imp_brain/conn_regions_reordered.txt'), 'r') as f:
  2304. final_order = [line.strip() for line in f.readlines()]
  2306. # Assign 'Subcortical' as network for subcortical regions (FIXED: Do this FIRST)
  2307. coords['network'] = coords.apply(
  2308. lambda row: row['network'] if row['type'] == 'cortical' else 'Subcortical',
  2309. axis=1
  2310. )
  2312. # Assign hemisphere for subcortical regions
  2313. coords['hemi'] = coords.apply(
  2314. lambda row: row['hemi'] if row['type'] == 'cortical' else (
  2315. row['name'].split('-')[-1] if row['name'].endswith(('-lh', '-rh')) else 'none'
  2316. ),
  2317. axis=1
  2318. )
  2320. # Align coordinates to final order
  2321. coords_aligned = (
  2322. coords.set_index('name')
  2323. .loc[final_order]
  2324. .reset_index()
  2325. )
  2327. # Convert to correlation matrices
  2328. corr_matrix_dwi, pval_matrix_dwi = long_to_square_matrix_dwi(corr_dwi, coords)
  2329. corr_matrix_rs, pval_matrix_rs = long_to_square_matrix_rs(corr_rs, coords)
  2331. # Reorder correlation matrices
  2332. corr_matrix_ordered_dwi = corr_matrix_dwi.loc[final_order, final_order]
  2333. corr_matrix_ordered_rs = corr_matrix_rs.loc[final_order, final_order]
  2335. # Set up network colors and names
  2336. network_colors = {
  2337. 'Vis': '#ff0000',
  2338. 'SomMot': '#00fa9a',
  2339. 'DorsAttn': '#0076BBFF',
  2340. 'SalVentAttn': '#00ced1',
  2341. 'Limbic': '#ffff00',
  2342. 'Cont': '#BBBBBB',
  2343. 'Default': '#ff69b4',
  2344. 'Subcortical': '#5AAE61FF'
  2345. }
  2347. structure_names = {
  2348. 'Vis': 'Vis',
  2349. 'SomMot': 'SomMot',
  2350. 'DorsAttn': 'DAN',
  2351. 'SalVentAttn': 'VAN',
  2352. 'Limbic': 'Limb',
  2353. 'Cont': 'Cont',
  2354. 'Default': 'DMN',
  2355. 'Subcortical': 'Subcort'
  2356. }
  2358. structure_names_full = {
  2359. 'Vis': 'Visual',
  2360. 'SomMot': 'Somatomotor',
  2361. 'DorsAttn': 'Dorsal Attention',
  2362. 'SalVentAttn': 'Salience/Ventral Attention',
  2363. 'Limbic': 'Limbic',
  2364. 'Cont': 'Control',
  2365. 'Default': 'Default Mode',
  2366. 'Subcortical': 'Subcortical'
  2367. }
  2369. # %%
  2370. # Plot heatmap: vertical
  2371. n = len(corr_matrix_ordered_dwi)
  2373. # Create figure with two subplots side by side
  2374. fig, axes = plt.subplots(2, 1, figsize=(16, 24), gridspec_kw={'height_ratios': [1, 1], 'hspace': 0.2})
  2376. cmap = sns.color_palette("seismic", as_cmap=True)
  2377. vmin, vmax = -0.45, 0.45
  2379. # Create lower triangle mask (hides lower triangle, shows upper triangle)
  2380. #upper_triangle_mask = np.tril(np.ones_like(corr_matrix_ordered_dwi, dtype=bool), k=-1)
  2381. lower_triangle_mask = np.triu(np.ones_like(corr_matrix_ordered_dwi, dtype=bool), k=1)
  2383. # Top
  2384. ax1 = axes[0]
  2386. rs_data = corr_matrix_ordered_rs.values.copy()
  2387. np.fill_diagonal(rs_data, 0)
  2389. sns.heatmap(
  2390. rs_data,
  2391. mask=lower_triangle_mask,
  2392. cmap=cmap,
  2393. center=0,
  2394. vmin=vmin, vmax=vmax,
  2395. square=True,
  2396. xticklabels=False,
  2397. yticklabels=False,
  2398. linewidths=0,
  2399. cbar=False,
  2400. ax=ax1,
  2401. rasterized=True
  2402. )
  2404. # Diagonal line for left plot (bottom-left to top-right)
  2405. ax1.plot([0, n], [0, n], color='black', linewidth=1, linestyle='-', alpha=0.7)
  2407. # Add rsMRI label
  2408. #ax1.text(-0.1, 0.5, 'rsMRI', transform=ax1.transAxes, fontsize=30, va='center', ha='right', rotation=90)
  2410. # Bottom
  2411. ax2 = axes[1]
  2413. dwi_data = corr_matrix_ordered_dwi.values.copy()
  2414. np.fill_diagonal(dwi_data, 0)
  2416. sns.heatmap(
  2417. dwi_data,
  2418. mask=lower_triangle_mask, # mask shows upper triangle
  2419. cmap=cmap,
  2420. center=0,
  2421. vmin=vmin, vmax=vmax,
  2422. square=True,
  2423. xticklabels=False,
  2424. yticklabels=False,
  2425. linewidths=0,
  2426. cbar=False,
  2427. ax=ax2,
  2428. rasterized=True
  2429. )
  2431. # Diagonal line for bottom plot
  2432. ax2.plot([0, n], [0, n], color='black', linewidth=1, linestyle='-', alpha=0.7)
  2434. # Add dwMRI label
  2435. #ax2.text(1.1, 0.5, 'dwMRI', transform=ax2.transAxes, fontsize=30, va='center', ha='left', rotation=270)
  2437. # Add shared colorbar
  2438. from matplotlib.cm import ScalarMappable
  2439. from matplotlib.colors import Normalize
  2441. sm = ScalarMappable(cmap=cmap, norm=Normalize(vmin=vmin, vmax=vmax))
  2442. sm.set_array([])
  2444. cbar_ax = fig.add_axes([0.3, 0.08, 0.4, 0.02])
  2445. cbar = fig.colorbar(sm, cax=cbar_ax, orientation='horizontal')
  2446. cbar.set_label('Pearson $r$', fontsize=50)
  2447. cbar.ax.tick_params(labelsize=30)
  2449. # Annotations
  2450. bottom_y = n + 1.5
  2451. left_x = -8
  2452. left_x -= 12
  2453. bar_height = 25
  2454. bar_width = 20
  2456. # Get network boundaries
  2457. network_boundaries = []
  2458. current_network = None
  2459. for i, region in enumerate(corr_matrix_ordered_dwi.index):
  2460. network = coords_aligned.loc[coords_aligned['name'] == region, 'network'].values[0]
  2461. if network != current_network:
  2462. if current_network is not None:
  2463. network_boundaries.append((start_idx, i, current_network))
  2464. start_idx = i
  2465. current_network = network
  2466. network_boundaries.append((start_idx, n, current_network))
  2468. # Apply to bottom plot
  2469. for start, end, network in network_boundaries:
  2470. color = network_colors[network]
  2472. # Bottom color bar
  2473. ax1.add_patch(plt.Rectangle(
  2474. (start, bottom_y - 2),
  2475. end - start,
  2476. bar_height,
  2477. facecolor=color,
  2478. clip_on=False
  2479. ))
  2481. # Left color bar
  2482. ax1.add_patch(plt.Rectangle(
  2483. (left_x, start),
  2484. bar_width,
  2485. end - start,
  2486. facecolor=color,
  2487. clip_on=False
  2488. ))
  2490. # Structure labels
  2491. mid_point = (start + end) / 2
  2492. # Left
  2493. #ax1.text(left_x + 17, mid_point - 2.1, structure_names.get(network, network), ha='right', va='center', fontsize=15, rotation=90)
  2494. # Bottom
  2495. #ax1.text(mid_point, bottom_y + bar_height/1.5 + 2, structure_names.get(network, network), ha='center', va='bottom', fontsize=15)
  2497. # Apply to dwMRI
  2498. for start, end, network in network_boundaries:
  2499. color = network_colors[network]
  2501. # Bottom color bar
  2502. ax2.add_patch(plt.Rectangle(
  2503. (start, bottom_y - 2),
  2504. end - start,
  2505. bar_height,
  2506. facecolor=color,
  2507. clip_on=False
  2508. ))
  2509. # Left color bar
  2510. ax2.add_patch(plt.Rectangle(
  2511. (left_x, start),
  2512. bar_width,
  2513. end - start,
  2514. facecolor=color,
  2515. clip_on=False
  2516. ))
  2518. # Structure labels for bottom
  2519. mid_point = (start + end) / 2
  2520. # Left
  2521. #ax2.text(left_x + 17, mid_point - 2.1, structure_names.get(network, network), ha='right', va='center', fontsize=15, rotation=90)
  2522. # Bottom
  2523. #ax2.text(mid_point, bottom_y + bar_height/1.5 + 2, structure_names.get(network, network), ha='center', va='bottom', fontsize=15)
  2525. # Grid lines and limits for LEFT plot
  2526. top_y = -10
  2527. for boundary in [b[0] for b in network_boundaries[1:]]:
  2528. # For upper triangle orientation
  2529. ax1.plot([boundary, boundary], [boundary, n], color='black', linewidth=0.4, alpha=1) # Vertical
  2530. ax1.plot([0, boundary], [boundary, boundary], color='black', linewidth=0.4, alpha=1) # Horizontal
  2532. ax1.set_xlim(left_x - 1, n + 0.8)
  2533. ax1.set_ylim(n + bar_height + 2, top_y - 2)
  2535. # Only draw left and bottom borders
  2536. ax1.plot([0, 0], [0, n], color='black', linewidth=1) # Left border
  2537. ax1.plot([0, n], [n, n], color='black', linewidth=1) # Bottom border
  2539. # Grid lines and limits for Rbottomplot
  2540. for boundary in [b[0] for b in network_boundaries[1:]]:
  2541. ax2.plot([boundary, boundary], [boundary, n], color='black', linewidth=0.4, alpha=1) # Vertical
  2542. ax2.plot([0, boundary], [boundary, boundary], color='black', linewidth=0.4, alpha=1) # Horizontal
  2544. ax2.set_xlim(-1, n + 0.8)
  2545. ax2.set_ylim(n + bar_height + 2, top_y - 2)
  2547. # Only draw left and bottom borders
  2548. ax2.plot([0, 0], [0, n], color='black', linewidth=1) # Left border
  2549. ax2.plot([0, n], [n, n], color='black', linewidth=1) # Bottom border
  2551. # Individual titles for each subplot
  2552. #ax1.set_title("rsMRI Feature Importance:\nSchaefer7n500p-IV Full Correlation", fontsize=36, pad=20)
  2553. #ax2.set_title("dwMRI Feature Importance:\nSchaefer7n500p-IV Streamline Count", fontsize=36, pad=20)
  2555. # Clean up all spines
  2556. for ax in [ax1, ax2]:
  2557. ax.spines['top'].set_visible(False)
  2558. ax.spines['right'].set_visible(False)
  2559. ax.spines['left'].set_visible(False)
  2560. ax.spines['bottom'].set_visible(False)
  2562. # Extract unique networks
  2563. unique_networks = sorted(set(
  2564. coords_aligned.loc[coords_aligned['name'].isin(corr_matrix_ordered_dwi.index), 'network']
  2565. ))
  2567. # Create legend elements
  2568. legend_elements = []
  2569. # Set up network order matching patches
  2570. unique_networks = [network for _, _, network in network_boundaries]
  2571. for network in unique_networks:
  2572. color = network_colors.get(network, 'gray')
  2573. label = structure_names_full.get(network, network)
  2574. legend_elements.append(
  2575. plt.Rectangle((0, 0), 1, 1, facecolor=color, edgecolor='black',
  2576. linewidth=0.5, label=label)
  2577. )
  2579. # Add legend to figure
  2580. legend = fig.legend(
  2581. handles=legend_elements,
  2582. loc='lower center',
  2583. ncol= 1, #min(len(unique_networks), 9), # Adjust columns
  2584. fontsize=40,
  2585. frameon=True,
  2586. fancybox=True,
  2587. framealpha=0.9,
  2588. bbox_to_anchor=(1.0, 0.6), #vertical, on the right side
  2589. #bbox_to_anchor=(0.5, -0.01), horizontal, one row
  2590. handletextpad=0.5,
  2591. columnspacing=1.0,
  2592. #title='Brain Networks',
  2593. title_fontsize=45
  2594. )
  2596. # Make room for the legend
  2597. plt.subplots_adjust(bottom=0.15)
  2599. # Save
  2600. save_fig(['png'], os.path.join(fig_path, 'Fig4c_connectomes_heatmaps_separate'))
  2601. plt.show()
  2603. # %%
  2604. '''
  2605. (x=0,y=0) -----------------(x=n,y=0)
  2606. | |
  2607. | |
  2608. | |
  2609. | |
  2610. (x=0,y=n) -----------------(x=n,y=n)
  2611. '''
  2612. '''
  2613. (0,0) → x=0, y=0 → the top-left corner of the matrix.
  2614. (n,0) → x=n, y=0 → the top-right corner of the matrix.
  2615. (0,n) → x=0, y=n → the bottom-left corner.
  2616. (n,n) → x=n, y=n → the bottom-right corner.
  2617. '''
  2619. # %% [markdown]
  2620. # # Figure 4d
  2622. # %% [markdown]
  2623. # sMRI barplots
  2625. # %%
  2626. # Filter to only show correlations > 0.2 in absolute value and p-value < 0.05
  2627. corr_smri = pd.read_csv(os.path.join(base_path, 'feature_imp/feature_imp_brain/struct_fs_aseg_volume_corr_with_g_stack.csv'))
  2628. plot_df_filtered = corr_smri[(corr_smri['Pearson r'].abs() > 0.195) & (corr_smri['p-value'] < 0.05)]
  2630. # Sort the filtered data
  2631. ordered = (plot_df_filtered
  2632. .assign(abs_corr = plot_df_filtered['Pearson r'].abs())
  2633. .sort_values('abs_corr', ascending=False)
  2634. ).sort_values('Pearson r', ascending=False)
  2636. # Create the figure with larger dimensions
  2637. n_bars = len(ordered)
  2638. plt.figure(figsize=(20, max(20, n_bars * 0.4)), dpi=300)
  2640. # Create color list - simple positive\negative coloring
  2641. colors = ['#D79C9CFF' if x > 0 else '#5F93ACFF' for x in ordered['Pearson r']]
  2643. # Create bars with slightly thicker bars
  2644. bars = plt.barh(ordered['Feature'].str.replace(' hemisphere', ''),
  2645. ordered['Pearson r'],
  2646. color=colors,
  2647. height=0.8)
  2649. # Add correlation values with larger font
  2650. for bar, corr in zip(bars, ordered['Pearson r']):
  2651. width = bar.get_width()
  2652. label_x = width - 0.01 if width > 0 else width + 0.02
  2653. plt.text(label_x, bar.get_y() + bar.get_height()/2,
  2654. f'{width:.2f}',
  2655. va='center', ha='right' if width > 0 else 'left',
  2656. color='white', fontsize=30)
  2658. # Increase all font sizes
  2659. #plt.title(f'sMRI Feature Importance:\nASEG Subcortical Volumetric Segmentation', pad=18, fontsize=45, x=-0.5)
  2660. plt.xlabel('Pearson $r$', labelpad=10, fontsize=45)
  2661. plt.xticks(fontsize=35)
  2662. plt.ylabel('')
  2663. plt.yticks(fontsize=37)
  2664. plt.gca().invert_yaxis()
  2666. # Adjust layout with more padding
  2667. plt.tight_layout(pad=3.0)
  2669. # Save with higher quality
  2670. plt.savefig(os.path.join(fig_path, 'Fig4d_fi-smri-BARS.svg'), dpi=350,bbox_inches='tight')
  2671. plt.show()
  2673. # %% [markdown]
  2674. # # Figure 5
  2676. # %% [markdown]
  2677. # g-body commonality analysis
  2679. # %%
  2680. # Prepare data for plotting
  2681. commonality_path = '/UK_BB/brainbody/commonality'
  2682. df = pd.read_csv(os.path.join(commonality_path, 'commonality_results_decompos_brain_vs_individual_body.csv'))
  2683. print(df.columns.to_list(), df['Brain MRI modality'].unique(), df['Body modality'].unique())
  2685. # Dataset 1: For Venn diagrams (stacked body modalities)
  2686. df_venn = pd.read_excel(os.path.join(commonality_path, 'commonality_results_decompos_brain_vs_body_stacked.xlsx'))
  2688. # Dataset 2: For stacked bar plots (individual body modalities)
  2689. df_bars = pd.read_csv(os.path.join(commonality_path, 'commonality_results_decompos_brain_vs_individual_body.csv'))
  2691. # Check brain modality names in both DataFrames
  2692. print("Brain modalities in df_venn (Venn diagram data):")
  2693. print(df_venn['Brain MRI modality'].unique())
  2694. print("\nBrain modalities in df_bars (bar plot data):")
  2695. print(df_bars['Brain MRI modality'].unique())
  2697. # Filter for Body physiology and composition stacked
  2698. body_modality = 'Body physiology and composition stacked'
  2699. body_df = df_venn[df_venn['Body modality'] == body_modality]
  2701. # Get all brain modalities
  2702. brain_modalities_venn = body_df['Brain MRI modality'].unique()
  2703. n_modalities = len(brain_modalities_venn)
  2705. # Prepare venn2 data
  2706. data_venn2 = []
  2707. color_palette = {
  2708. 'brain': '#4A7169FF',#'#4A7169FF'
  2709. 'body': '#A8554EFF' #'#A8554EFF'
  2710. }
  2711. for brain_mod in brain_modalities_venn:
  2712. row = body_df[body_df['Brain MRI modality'] == brain_mod].iloc[0]
  2713. prop_body_exp_by_brain = round((row['Common variance'] / (row['Unique variance: body'] + row['Common variance'])) * 100, 2)
  2714. data_venn2.append((
  2715. #row['Perc unique brain'],
  2716. #row['Perc unique body'],
  2717. #row['Perc common'],
  2718. row['Unique variance: brain'],
  2719. row['Unique variance: body'],
  2720. row['Common variance'],
  2721. color_palette['brain'],
  2722. color_palette['body'],
  2723. brain_mod.replace('Brain ', '').replace('3 brain', '3'),
  2724. prop_body_exp_by_brain
  2725. ))
  2727. # Create a mapping from Venn names to bar plot names
  2728. brain_mod_mapping = {
  2729. '3 brain MRI modalities stacked': '3 Brain MRI Modalities Stacked',
  2730. 'Brain dwMRI stacked': 'Brain dwMRI Stacked',
  2731. 'Brain rsMRI stacked': 'Brain rsMRI Stacked',
  2732. 'Brain sMRI stacked': 'Brain sMRI Stacked'
  2733. }
  2735. # Get the brain modality names for bar plots
  2736. brain_modalities_bars = [brain_mod_mapping.get(mod, mod) for mod in brain_modalities_venn]
  2738. print(f"\nVenn brain modalities: {brain_modalities_venn}")
  2739. print(f"Bar plot brain modalities: {brain_modalities_bars}")
  2741. # Get the sorted order from the first brain modality
  2742. first_brain_mod_bars = brain_modalities_bars[0]
  2743. df_first = df_bars[df_bars['Brain MRI modality'] == first_brain_mod_bars].copy()
  2744. print(f"\nFirst brain modality for sorting: '{first_brain_mod_bars}'")
  2745. print(f"Number of rows found: {len(df_first)}")
  2747. df_first.sort_values(by='Unique variance: brain', ascending=False, inplace=True)
  2748. sorted_body_order = df_first['Body modality'].tolist()
  2749. print(f"Sorted body order length: {len(sorted_body_order)}")
  2751. # Create a dictionary to store data for each subplot
  2752. plot_data = {}
  2753. for venn_mod, bars_mod in zip(brain_modalities_venn, brain_modalities_bars):
  2754. print(f"Processing: Venn='{venn_mod}' -> Bars='{bars_mod}'")
  2756. # Filter dataframe for the current brain modality
  2757. df_mod = df_bars[df_bars['Brain MRI modality'] == bars_mod].copy()
  2758. print(f" Rows found: {len(df_mod)}")
  2760. # Create a categorical type with the desired order
  2761. df_mod['Body modality'] = pd.Categorical(df_mod['Body modality'],
  2762. categories=sorted_body_order,
  2763. ordered=True)
  2764. # Sort by the predefined order
  2765. df_mod.sort_values('Body modality', inplace=True)
  2767. # Store the data for this subplot
  2768. plot_data[venn_mod] = {
  2769. 'body_names': df_mod['Body modality'].tolist(),
  2770. 'y1_unique_brain': df_mod['Unique variance: brain'].values,
  2771. 'y2_common': df_mod['Common variance'].values,
  2772. 'y3_unique_body': df_mod['Unique variance: body'].values,
  2773. # Store percentages for annotations
  2774. 'perc_u_brain': df_mod['Perc unique brain'].values,
  2775. 'perc_common': df_mod['Perc common'].values,
  2776. 'perc_u_body': df_mod['Perc unique body'].values,
  2777. 'prop_body_exp_by_brain': df_mod['g-body explained by brain'].values,
  2778. 'prop_body_exp_by_brain': df_mod['g-body explained by brain'].values
  2779. }
  2781. # %%
  2782. # Plot
  2783. fig = plt.figure(figsize=(6*n_modalities, 24))
  2785. # Define manual positions
  2786. positions = []
  2787. # Top row - Venn diagrams
  2788. for i in range(n_modalities):
  2789. left = 0.05 + i * (0.9/n_modalities) + i*0.04
  2790. positions.append([left, 0.61, 0.8/n_modalities, 0.25]) # [left, bottom, width, height]
  2792. # Bottom row - Bar plots
  2793. for i in range(n_modalities):
  2794. left = 0.05 + i * (0.9/n_modalities) + i*0.04
  2795. positions.append([left, 0.05, 0.8/n_modalities, 0.5])
  2797. # Create axes
  2798. axes = []
  2799. for pos in positions:
  2800. axes.append(fig.add_axes(pos))
  2802. # Define colors for stacked bars
  2803. color_unique_brain = '#4A7169FF'
  2804. color_common = '#c2ada6'
  2805. color_unique_body = '#A8554EFF'
  2807. # Plot Venn diagrams (top row - axes[0] to axes[n_modalities-1])
  2808. for i, (u_brain, u_body, c_common, color_brain, color_body, label, prop_body_exp) in enumerate(data_venn2):
  2809. ax = axes[i] # Use first n_modalities axes for Venn diagrams
  2811. # Create venn2 diagram
  2812. subset_sizes = (u_brain, u_body, c_common)
  2813. #subset_sizes = ( row['Unique variance: brain'], row['Unique variance: body'], row['Common variance'])
  2814. venn = venn2(subsets=subset_sizes, set_labels=('', ''),
  2815. set_colors=(color_brain, color_body), ax=ax)
  2817. ax.set_ylim(-1.2, 1.2)
  2818. ax.set_aspect(0.8) #'equal'
  2819. ax.axis('off')
  2821. # Apply styling
  2822. for patch in venn.patches:
  2823. if patch:
  2824. patch.set_alpha(0.7)
  2825. patch.set_linewidth(1.5)
  2826. patch.set_edgecolor("black")
  2828. # Format labels
  2829. #for text in venn.subset_labels:
  2830. for j, text in enumerate(venn.subset_labels):
  2831. if text:
  2832. val = float(text.get_text())
  2833. text.set_text(f'{val*100:.1f}')
  2834. text.set_fontsize(30)
  2836. # Adjust the leftmost label (index 0 corresponds to "10" region, i.e. left circle only)
  2837. if j == 0: # left-only subset
  2838. x, y = text.get_position()
  2839. text.set_position((x + 0.04, y)) # shift right by 0.05 units
  2841. # Adjust Venn diagram annotations
  2842. if i == 1:
  2843. if j == 0:
  2844. x, y = text.get_position()
  2845. text.set_position((x - 0.035, y))
  2847. if i == 3:
  2848. if j == 0:
  2849. x, y = text.get_position()
  2850. text.set_position((x - 0.035, y))
  2852. # Add title and annotation using figure coordinates
  2853. titles = [ "Composite Brain",
  2854. "dwMRI",
  2855. "rsMRI",
  2856. "sMRI" ]
  2858. # Title
  2859. title_x = positions[i][0] + positions[i][2]/2 # Center of the Venn axes
  2860. fig.text(title_x, 0.83, titles[i], fontsize=45, ha='center', va='top')
  2862. # Annotation
  2863. annotation_text = f"{prop_body_exp:.1f}%"
  2864. fig.text(title_x, 0.65, annotation_text, ha='center', va='center',
  2865. fontsize=35,
  2866. bbox=dict(boxstyle="round,pad=0.3", facecolor='white', alpha=0.7))
  2869. for idx, brain_mod in enumerate(brain_modalities_venn):
  2870. ax = axes[n_modalities + idx]
  2871. data = plot_data[brain_mod]
  2872. y_labels = data['body_names']
  2873. # Stacked barplots (bottom row - axes[n_modalities] to axes[2*n_modalities-1])
  2874. y_labels = [' '.join(word.upper() if word.upper() in ['MRI', 'DXA', 'ECG'] else word.lower() if i > 0 else word.capitalize() for i, word in enumerate(label.split()))
  2875. for label in y_labels]
  2876. rename_dict = {
  2877. 'Bone densitometry (of heel)': 'Bone densitometry (heel)',
  2878. 'Abdominal organ composition by MRI': 'Abdominal organ composition (MRI)',
  2879. 'Body composition by DXA': 'Body composition (DXA)',
  2880. 'Abdominal composition by MRI': 'Abdominal composition (MRI)',
  2881. 'Bone size, mineral and density by DXA': 'BMD (DXA)',
  2882. 'Body composition by impedance': 'Body Composition (impedance)'
  2883. }
  2885. y_labels = [rename_dict.get(lbl, lbl) for lbl in y_labels]
  2887. y1 = data['y1_unique_brain']*100
  2888. y2 = data['y2_common']*100
  2889. y3 = data['y3_unique_body']*100
  2890. prop_body_exp_by_brain = data['prop_body_exp_by_brain']
  2891. bar_height = 0.8
  2893. # Create stacked horizontal bars
  2894. ax.barh(y_labels, y1, height=bar_height, color=color_unique_brain, alpha=0.8, edgecolor='black', linewidth=1.5)
  2895. ax.barh(y_labels, y2, height=bar_height, left=y1, color=color_common, alpha=0.8, edgecolor='black', linewidth=1.5)
  2896. ax.barh(y_labels, y3, height=bar_height, left=y1+y2, color=color_unique_body, alpha=0.8, edgecolor='black', linewidth=1.5)
  2898. # Add percentage annotations on each segment
  2899. for i, (u_brain, common, u_body) in enumerate(zip(y1, y2, y3)):
  2900. if u_brain > 0.001:
  2901. ax.text(u_brain / 2, i, f'{u_brain:.1f}', ha='center', va='center',
  2902. color='white', fontsize=22)
  2903. if common > 0.001:
  2904. if common < 0.9: # Small common variance - to the left
  2905. # Position at the start of the common segment, shifted slightly right
  2906. ax.text(u_brain + common / 2, i, f'{common:.1f}', ha='right', va='center',
  2907. color='white', fontsize=22)
  2908. elif common < 2 and common > 1:
  2909. ax.text(u_brain + common / 2, i, f'{common:.1f}', ha='right', va='center', color='white', fontsize=20)
  2910. else: # Normal/large common variance - center
  2911. ax.text(u_brain + common / 2, i, f'{common:.1f}', ha='center', va='center', color='white', fontsize=20)
  2912. if u_body: # > 0.001:
  2913. #ax.text(u_brain + common + u_body / 2, i, f'{m_u_body:.2f}', ha='center', va='center', color='white', fontsize=20)
  2914. pass
  2915. prop_text = f'{prop_body_exp_by_brain[i]:.1f}%'
  2916. ax.text(1.25, i, prop_text, ha='right', va='center',
  2917. fontsize=25, transform=ax.get_yaxis_transform(), color='black', bbox=dict(boxstyle="round,pad=0.2", facecolor='white', alpha=0.7))
  2919. # Format the subplot
  2920. ax.tick_params(axis='y', labelsize=30, length=0)
  2921. ax.tick_params(axis='x', labelsize=45)
  2922. ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%.1f'))
  2923. max_total = max(y1 + y2 + y3)
  2924. ax.set_xlim(0, max_total)
  2925. ax.set_xticks([0, max_total])
  2927. # Put ticks at 0, half, and max
  2928. import matplotlib.ticker as mticker
  2929. half_point = max_total / 2
  2930. ax.xaxis.set_major_locator(mticker.FixedLocator([0, half_point, max_total]))
  2932. if idx > 0:
  2933. ax.set_yticklabels([])
  2935. for spine in ['top', 'right', 'left']:
  2936. ax.spines[spine].set_visible(False)
  2938. # Add common x-axis label for bar plots
  2939. fig.text(0.53, -0.03, 'Variance of Cognition Explained (%)', ha='center', fontsize=50)
  2940. fig.suptitle('Commonality Analysis\nbetween Body Phenotypes and Brain Markers\nin Explaining Cognition (%)', fontsize=70, y=1)
  2943. def make_circle_legend(legend, orig_handle, xdescent, ydescent, width, height, fontsize):
  2944. return plt.Circle((width/2, height/2), min(width, height))
  2946. legend_elements = [
  2947. plt.Circle((0,0), 0.5, facecolor=color_unique_brain, alpha=0.9, label='Unique to Brain'),
  2948. plt.Circle((0,0), 0.5, facecolor=color_common, alpha=0.9, label='Common Variance'),
  2949. plt.Circle((0,0), 0.5, facecolor=color_unique_body, alpha=0.9, label='Unique to Body'),
  2950. plt.Rectangle((0,0), 1, 1, facecolor='none', edgecolor='black', lw=2,
  2951. label='Body-related cognitive variance\nexplained by each brain marker (%)')]
  2953. leg = fig.legend(handles=legend_elements, loc='lower center', ncol=4,
  2954. fontsize=38, frameon=True, #bbox_to_anchor=(0.5, -0.005),
  2955. fancybox=True,
  2956. bbox_to_anchor=(0.5, 0.54),
  2957. handler_map={plt.Circle: HandlerPatch(patch_func=make_circle_legend)},
  2958. handletextpad=1.0, handleheight=1.0,
  2959. columnspacing=1.5, labelspacing=1.5,
  2960. borderpad=0.5
  2961. )
  2963. # Control the legend frame linewidth and color
  2964. frame = leg.get_frame()
  2965. frame.set_linewidth(2)
  2966. frame.set_edgecolor("black"),
  2967. frame.set_facecolor(((216/255, 216/255, 216/255, 1)))
  2969. # Save
  2970. save_paths = [
  2971. os.path.join(fig_path, 'final/Fig5'),
  2972. os.path.join(fig_path, 'Fig5_ca_g_body_exp_by_brain_stacked_bars_plus_venn_marginalR2'),]
  2973. for path in save_paths:
  2974. save_fig(['png', 'svg'], path)
  2975. plt.show()
  2977. # %%
  2978. # Venn 2 for Fig1
  2979. fig, ax = plt.subplots(figsize=(8, 8))
  2981. u_brain, u_body, c_common, color_brain, color_body, label, prop_body_exp = data_venn2[0]
  2983. subset_sizes = (u_body, u_brain, c_common)
  2984. venn = venn2(subsets=subset_sizes, set_labels=('Body','Brain'), # add labels
  2985. set_colors=(color_body, '#0D5581FF'), ax=ax)
  2987. ax.set_ylim(-1.2, 1.2)
  2988. ax.set_aspect(0.9)
  2989. ax.axis('off')
  2991. # Apply styling
  2992. for patch in venn.patches:
  2993. if patch:
  2994. patch.set_alpha(0.7)
  2995. patch.set_linewidth(1.5)
  2996. patch.set_edgecolor("black")
  2998. # Remove text inside circles
  2999. for text in venn.subset_labels:
  3000. if text:
  3001. text.set_text('')
  3003. # Labels
  3004. for i, text in enumerate(venn.set_labels):
  3005. if text:
  3006. text.set_fontsize(32)
  3007. label_colors = [color_body, '#0D5581FF']
  3008. text.set_color(label_colors[i])
  3010. # Position adjustments
  3011. if i == 0: # "Brain" label - move left
  3012. x, y = text.get_position()
  3013. text.set_position((x - 0.1, y + 1.2))
  3014. elif i == 1: # "Body" label - move right
  3015. x, y = text.get_position()
  3016. text.set_position((x + 0.2, y + 1.3))
  3017. # Save
  3018. save_fig(['svg'], os.path.join(fig_path, 'Fig1d_venn2'))
  3019. plt.show()
  3021. # %% [markdown]
  3022. # # Figure 6
  3024. # %% [markdown]
  3025. # ## Figure 6a
  3027. # %% [markdown]
  3028. # g-age commonality analysis: venn3 diagrams
  3030. # %%
  3031. # Set the path
  3032. commonality_path = '/UK_BB/brainbody/commonality'
  3034. # %%
  3035. # Prepare data for plotting
  3036. commonality_path = 'Volumes/sci-psy-narun/IBu/UK_BB/brainbody/commonality'
  3038. df = pd.read_csv(
  3039. os.path.join(commonality_path, 'commonality_results_body_plus_brain_age_3v.csv')
  3040. )
  3042. # Extract raw variance components
  3043. u_body = df['Unique: Body'].iloc[0]
  3044. u_brain = df['Unique: Brain'].iloc[0]
  3045. u_age = df['Unique: Age'].iloc[0]
  3047. c_body_brain = df['Common: Body & Brain'].iloc[0]
  3048. c_body_age = df['Common: Body & Age'].iloc[0]
  3049. c_brain_age = df['Common: Brain & Age'].iloc[0]
  3051. c_all_three = df['Common: All Three'].iloc[0]
  3053. # Pack into Venn3 order
  3054. subset_values = (
  3055. u_body, # 100
  3056. u_brain, # 010
  3057. c_body_brain, # 110
  3058. u_age, # 001
  3059. c_body_age, # 101
  3060. c_brain_age, # 011
  3061. c_all_three # 111
  3062. )
  3064. # Total R2
  3065. total_r2 = sum(subset_values)
  3066. print("Total R2:", round(total_r2, 4))
  3068. # Convert to percentages for plotting (Venn3 geometry)
  3069. subset_percent = tuple(v / total_r2 * 100 for v in subset_values)
  3070. print("%:", subset_percent)
  3072. # %%
  3073. # Venn 3
  3074. fig, ax = plt.subplots(figsize=(10, 10))
  3076. body = "#9E2F2584" #A8554E
  3077. brain = "#0076C0FF" #4A7169
  3078. age = "#2B811AC0" #C2ADA6
  3080. body_label = "#74251EFF" #A8554E
  3081. brain_label = "#084972FF" #4A7169
  3082. age_label = "#104905C0" #C2ADA6
  3084. # Plot using percentages
  3085. venn = venn3(
  3086. subsets=subset_percent,
  3087. set_labels=('Body', 'Brain', 'Age'),
  3088. set_colors=(body, brain, age),
  3089. ax=ax
  3090. )
  3092. ax.set_aspect(0.9)
  3094. for patch in venn.patches:
  3095. if patch:
  3096. patch.set_alpha(0.6)
  3097. patch.set_edgecolor('black')
  3098. patch.set_linewidth(1.5)
  3100. # Replace percentage labels with raw variance values
  3101. for text, real_value in zip(venn.subset_labels, subset_values):
  3102. if text:
  3103. text.set_text(f"{real_value*100:.1f}")
  3104. text.set_fontsize(25)
  3105. text.set_color('white')
  3107. # Set labels ("Body", "Brain", "Age")
  3108. for i, text in enumerate(venn.set_labels):
  3109. if text:
  3110. text.set_fontsize(40)
  3111. #text.set_fontweight('bold')
  3112. label_colors = [body_label, brain_label, age_label]
  3113. text.set_color(label_colors[i]) # 'black'
  3115. # Get label positions
  3116. for region in ['100','010','001','110','101','011','111']:
  3117. lab = venn.get_label_by_id(region)
  3118. if lab:
  3119. print(region, lab.get_position())
  3121. # Move Body-only label
  3122. lab = venn.get_label_by_id('100')
  3123. if lab:
  3124. x, y = lab.get_position()
  3125. lab.set_position((x + 0.01, y + 0.02))
  3127. # Move Brain-only label
  3128. #lab = venn.get_label_by_id('010')
  3129. #if lab:
  3130. #x, y = lab.get_position()
  3131. #lab.set_position((x + 0.04, y))
  3133. # Move Age-only label
  3134. #lab = venn.get_label_by_id('001')
  3135. #if lab:
  3136. # x, y = lab.get_position()
  3137. #lab.set_position((x, y + 0.03))
  3139. # Move Body ∩ Brain
  3140. lab = venn.get_label_by_id('110')
  3141. if lab:
  3142. lab.set_position((-0.0197, 0.4))
  3144. # Move Body ∩ Age
  3145. lab = venn.get_label_by_id('101')
  3146. if lab:
  3147. x, y = lab.get_position()
  3148. lab.set_position((x - 0.03, y + 0.02))
  3150. # Move Brain ∩ Age
  3151. lab = venn.get_label_by_id('011')
  3152. if lab:
  3153. x, y = lab.get_position()
  3154. lab.set_position((x - 0.005, y - 0.05))
  3156. # Move triple overlap
  3157. #lab = venn.get_label_by_id('111')
  3158. #if lab:
  3159. #x, y = lab.get_position()
  3160. #lab.set_position((x - 0.03, y + 0.02))
  3162. # Add text on the right side of Venn diagrams
  3163. text_x_position = 1.0
  3164. text_y_start = 0.7
  3166. # Define items as (label, color, value)
  3167. # value=None means no percentage on that line
  3168. text_items = [
  3169. ('Age-related cognitive variance', None, None, None),
  3170. ('explained by:', None, None, None),
  3172. # Add a blank line marker to create extra spacing
  3173. ('', None, None, None),
  3175. ('Composite body:', '#74251EFF', '81.4%', 0.48),
  3176. ('Composite brain:', '#084972FF', '87.2%', 0.49),
  3177. ('Composite body or brain:', "#564464", '96.8%', 0.73),
  3178. ('Composite body and brain:', "#65746a", '71.7%', 0.77),
  3179. ]
  3180. # Vertical spacing
  3181. header_spacing = 0.06
  3182. line_spacing = 0.09
  3183. extra_gap_after_header = -0.05
  3185. y_pos = text_y_start
  3187. for i, (label, color, value, x_offset) in enumerate(text_items):
  3189. # Insert extra gap after the header block
  3190. if i == 3: # after the blank line
  3191. y_pos -= extra_gap_after_header
  3193. # Skip drawing empty lines
  3194. if label == '' and value is None:
  3195. y_pos -= line_spacing
  3196. continue
  3198. # Draw the label
  3199. if color:
  3200. fig.text(text_x_position, y_pos, label,
  3201. ha='left', va='center', fontsize=40,
  3202. bbox=dict(boxstyle="round,pad=0.2",
  3203. facecolor='white',
  3204. edgecolor=color,
  3205. linewidth=2.5),
  3206. transform=fig.transFigure)
  3207. else:
  3208. fig.text(text_x_position, y_pos, label,
  3209. ha='left', va='center', fontsize=40,
  3210. transform=fig.transFigure)
  3212. # Draw the percentage immediately after the label
  3213. if value and x_offset:
  3214. fig.text(text_x_position + x_offset, y_pos, # adjust horizontal offset
  3215. value,
  3216. ha='left', va='center', fontsize=40,
  3217. transform=fig.transFigure)
  3219. y_pos -= (header_spacing if i < 3 else line_spacing)
  3222. plt.title('Commonality Analysis\namong Age and Composite Body and Brain Markers\nin Explaining Cognition (%)', fontsize=55, y=1.2, x=0.9)
  3223. save_fig(['png', 'svg'], os.path.join(fig_path, 'Fig6b_3v_commonality_body_plus_brain_plus_age'))
  3224. plt.show()
  3226. # %%
  3227. # Venn 3
  3228. fig, ax = plt.subplots(figsize=(8, 8))
  3230. body = "#9E2F2584"
  3231. brain = "#0D5581FF"
  3232. age = "#2B811AC0"
  3234. body_label = "#74251EFF"
  3235. brain_label = "#084972FF"
  3236. age_label = "#104905C0"
  3238. # Plot Venn diagram
  3239. venn = venn3(
  3240. subsets=subset_percent,
  3241. set_labels=('Body', 'Brain', 'Age'),
  3242. set_colors=(body, brain, age),
  3243. ax=ax
  3244. )
  3246. # Remove axes completely
  3247. ax.axis('off')
  3249. # Remove all text labels (percentages inside circles)
  3250. for text in venn.subset_labels:
  3251. if text:
  3252. text.set_text('')
  3254. # Keep only the circle labels but make them smaller
  3255. for i, text in enumerate(venn.set_labels):
  3256. if text:
  3257. text.set_fontsize(40)
  3258. label_colors = [body_label, brain_label, age_label]
  3259. text.set_color(label_colors[i])
  3261. # Adjust circle transparency and edges
  3262. for patch in venn.patches:
  3263. if patch:
  3264. patch.set_alpha(0.6)
  3265. patch.set_edgecolor('black')
  3266. patch.set_linewidth(1.5)
  3268. # Save minimal figure
  3269. plt.savefig(os.path.join(fig_path, 'fig1-prep/Fig1d_venn3.svg'),
  3270. bbox_inches="tight",
  3271. pad_inches=0.5,
  3272. transparent=False,
  3273. facecolor="w",
  3274. edgecolor='w')
  3276. plt.show()
  3278. # %% [markdown]
  3279. # ## Figure 6b
  3281. # %% [markdown]
  3282. # g-age scatterplots
  3284. # %%
  3285. # Plot scatterplots for age ~ g
  3286. body_age = pd.read_csv(os.path.join(commonality_path, 'g_obs_pred_body_with_age.csv'))
  3287. brain_age = pd.read_csv(os.path.join(commonality_path, 'g_obs_pred_allmri_with_age.csv'))
  3288. scatter_kwargs = dict(
  3289. s=50,
  3290. alpha=0.9,
  3291. edgecolors='black',
  3292. linewidth=0.2
  3293. )
  3295. fig, axs = plt.subplots(1, 3, figsize=(35, 12))
  3297. # ---------------------------------------------------------------------
  3298. # Age vs observed g
  3299. x = pd.to_numeric(body_age['Age'], errors='coerce').astype(float)
  3300. y = pd.to_numeric(body_age['g_obs_test'], errors='coerce').astype(float)
  3301. axs[0].scatter(x, y, color="#ac87a0", **scatter_kwargs)
  3303. mask = np.isfinite(x) & np.isfinite(y)
  3304. if mask.sum() >= 2:
  3305. z = np.polyfit(x[mask], y[mask], 1)
  3306. p = np.poly1d(z)
  3307. xs = np.linspace(x[mask].min(), x[mask].max(), 200)
  3308. axs[0].plot(xs, p(xs), alpha=0.8, linewidth=2, color="#BE4A47FF")
  3310. axs[0].set_xlim(x.min(), x.max()+5)
  3311. y_min = y.min()
  3312. y_max = y.max()
  3313. tick_interval = 1
  3314. y_max_rounded = np.ceil(y_max / tick_interval) * tick_interval
  3315. axs[0].set_ylim(y_min, y_max_rounded)
  3316. ###
  3317. axs[0].set_ylim(y.min()-1, y.max()+1)
  3318. axs[0].set_title(r'$ĝ_{/text{observed}}$ ~ Age', fontsize=65, y=1.1, pad=15)#, bbox=dict(facecolor='white', alpha=0.9, edgecolor="#C98E8CFF", linewidth=2, boxstyle='round,pad=0.3'))
  3319. axs[0].set_xlabel('Age (years)', fontsize=55)
  3320. axs[0].set_ylabel('$g$-factor\nderived from ESEM ($z$)', fontsize=45)
  3321. axs[0].tick_params(axis='x', labelsize=42)
  3322. axs[0].tick_params(axis='y', labelsize=42)
  3323. axs[0].yaxis.set_major_locator(MultipleLocator(1))
  3324. # -------------------------------------------------------------------------
  3325. # Age vs predicted g (Body)
  3326. x = pd.to_numeric(body_age['Age'], errors='coerce').astype(float)
  3327. y = pd.to_numeric(body_age['g_pred_body_test'], errors='coerce').astype(float)
  3328. axs[1].scatter(x, y, color="#C98E8CFF", **scatter_kwargs)
  3330. mask = np.isfinite(x) & np.isfinite(y)
  3331. if mask.sum() >= 2:
  3332. z = np.polyfit(x[mask], y[mask], 1)
  3333. p = np.poly1d(z)
  3334. xs = np.linspace(x[mask].min(), x[mask].max(), 200)
  3335. axs[1].plot(xs, p(xs), alpha=0.8, linewidth=2, color="#BE4A47FF")
  3337. axs[1].set_xlim(x.min(), x.max()+5)
  3338. y_min = y.min()
  3339. y_max = y.max()
  3340. tick_interval = 1
  3341. y_max_rounded = np.ceil(y_max / tick_interval) * tick_interval
  3342. axs[1].set_ylim(y_min, y_max_rounded)
  3343. ###
  3344. axs[1].set_ylim(y.min()-1, y.max()+1)
  3345. axs[1].set_title(r'$ĝ_{/text{body}}$ ~ Age', fontsize=65, y=1.1, pad=15)#, bbox=dict(facecolor='white', alpha=0.9, edgecolor="#0d648f", linewidth=2, boxstyle='round,pad=0.3'))
  3346. axs[1].set_xlabel('Age (years)', fontsize=55)
  3347. axs[1].set_ylabel('$g$-factor\npredicted from body ($z$)', fontsize=45)
  3348. axs[1].tick_params(axis='x', labelsize=42)
  3349. axs[1].tick_params(axis='y', labelsize=42)
  3350. axs[1].yaxis.set_major_locator(MultipleLocator(1))
  3351. # ----------------------------------------------------------------------
  3352. # Age vs predicted cognition (Brain)
  3353. x = pd.to_numeric(brain_age['Age'], errors='coerce').astype(float)
  3354. y = pd.to_numeric(brain_age['g_pred_allmri_test'], errors='coerce').astype(float)
  3356. axs[2].scatter(x, y, color="#6A659999", **scatter_kwargs)
  3357. mask = np.isfinite(x) & np.isfinite(y)
  3358. if mask.sum() >= 2:
  3359. z = np.polyfit(x[mask], y[mask], 1)
  3360. p = np.poly1d(z)
  3361. xs = np.linspace(x[mask].min(), x[mask].max(), 200)
  3362. axs[2].plot(xs, p(xs), alpha=0.8, linewidth=2, color="#BE4A47FF")
  3364. axs[2].set_xlim(x.min(), x.max()+5)
  3365. y_min = y.min() - 1
  3366. y_max = y.max()
  3367. tick_interval = 1
  3368. y_max_rounded = np.ceil(y_max / tick_interval) * tick_interval
  3369. axs[2].set_ylim(y_min, y_max_rounded)
  3370. ###
  3371. axs[2].set_title(r'$ĝ_{/text{brain}}$ ~ Age', fontsize=65, y=1.1, pad=15)#, bbox=dict(facecolor='white', alpha=0.9, edgecolor="#f8b976", linewidth=2, boxstyle='round,pad=0.3'))
  3372. axs[2].set_xlabel('Age (years)', fontsize=55)
  3373. axs[2].set_ylabel('$g$-factor\npredicted from brain ($z$)', fontsize=50)
  3374. axs[2].tick_params(axis='x', labelsize=42)
  3375. axs[2].tick_params(axis='y', labelsize=42)
  3377. for ax in axs:
  3378. ax.spines['top'].set_visible(False)
  3379. ax.spines['right'].set_visible(False)
  3381. plt.tight_layout()
  3382. plt.subplots_adjust(wspace=0.34)
  3383. plt.savefig(os.path.join(fig_path, 'Fig6b_g_age_scatter.png'),
  3384. bbox_inches="tight",
  3385. pad_inches=1,
  3386. transparent=False,
  3387. facecolor="w",
  3388. edgecolor='w',
  3389. orientation='landscape',
  3390. format='png')
  3392. plt.show()
  3394. # %%
  3395. from scipy.stats import pearsonr, linregress
  3397. body_age = pd.read_csv(os.path.join(commonality_path, 'g_obs_pred_body_with_age.csv'))
  3398. brain_age = pd.read_csv(os.path.join(commonality_path, 'g_obs_pred_allmri_with_age.csv'))
  3400. def compute_metrics(x, y, name):
  3401. mask = ~(x.isna() | y.isna())
  3402. x_clean = x[mask]
  3403. y_clean = y[mask]
  3405. # Pearson r and p-value
  3406. r, p = pearsonr(x_clean, y_clean)
  3408. # R² from linear regression (same as r²)
  3409. slope, intercept, r_value, p_value, std_err = linregress(x_clean, y_clean)
  3411. print(f"\n{name}:")
  3412. print(f" r = {r:.3f}")
  3413. print(f" r² = {r**2:.3f}")
  3414. print(f" p = {p:.2e}")
  3416. # Compute for all three
  3417. compute_metrics(body_age['Age'], body_age['g_obs_test'], "Age vs observed g")
  3418. compute_metrics(body_age['Age'], body_age['g_pred_body_test'], "Age vs body-predicted g")
  3419. compute_metrics(brain_age['Age'], brain_age['g_pred_allmri_test'], "Age vs brain-predicted g")
  3421. # %% [markdown]
  3422. # # Figure 7
  3424. # %% [markdown]
  3425. # g-age commonality analysis: stacked bar plots
  3427. # %%
  3428. # Config
  3429. commonality_path = '/UK_BB/brainbody/commonality'
  3430. file_path = os.path.join(commonality_path,
  3431. 'commonality_results_decompos_brain_body_vs_age_sorted_stack_ind_combined_SORTED_BY_AGE_RENAMED.xlsx')
  3433. # Load data
  3434. df = pd.read_excel(file_path)
  3435. # Individual modalities (for bar plots)
  3436. df_individual = df[~df['Modality'].str.endswith('Stacked')].copy()
  3438. # Modality lists and order
  3439. modalities_to_plot = [
  3440. 'Brain dwMRI Stacked',
  3441. '3 Brain MRI Modalities Stacked',
  3442. 'Brain & Body Stacked',
  3443. 'Body Physiology Stacked',
  3444. 'Brain sMRI Stacked',
  3445. 'Brain rsMRI Stacked'
  3446. ]
  3448. modality_order = [
  3449. '3 Brain MRI Modalities Stacked',
  3450. 'Brain dwMRI Stacked',
  3451. 'Brain sMRI Stacked',
  3452. 'Brain & Body Stacked',
  3453. 'Brain rsMRI Stacked',
  3454. 'Body Physiology Stacked'
  3455. ]
  3457. # Colors
  3458. # For grouped bar plots
  3459. modality_colors = {
  3460. 'sMRI': '#f8b976',
  3461. 'dwMRI': '#0d648f',
  3462. 'rsMRI': '#4B6F5A99',
  3463. 'Body': '#C98E8CFF'
  3464. }
  3466. # Segment colors (shared)
  3467. color_unique_age = '#A8554EFF'
  3468. color_common = '#c2ada6'
  3469. color_unique_brain_body = '#4A7169FF'
  3471. # Case-insensitive mapping for individual modalities
  3472. lowercase_modality_mapping = {name.lower(): name for name in df_individual['Modality'].unique()}
  3474. modality_groups = {}
  3475. for group_name, mod_list in [('sMRI', modalities_smri),
  3476. ('dwMRI', modalities_dwmri),
  3477. ('rsMRI', modalities_rsmri),
  3478. ('Body', modalities_body)]:
  3479. modality_groups[group_name] = []
  3480. for mod in mod_list:
  3481. renamed_name = modality_names[mod]
  3482. lower_renamed = renamed_name.lower()
  3483. if lower_renamed in lowercase_modality_mapping:
  3484. modality_groups[group_name].append(lowercase_modality_mapping[lower_renamed])
  3485. else:
  3486. modality_groups[group_name].append(renamed_name)
  3488. print('dwMRI', len(modality_groups['dwMRI']))
  3489. print('rsMRI', len(modality_groups['rsMRI']))
  3490. print('sMRI', len(modality_groups['sMRI']))
  3491. print('Body', len(modality_groups['Body']))
  3493. # %% [markdown]
  3494. # ## Figure 7a
  3496. # %%
  3497. # Top Row - dwMRI and rsMRI
  3498. fig1, axs1 = plt.subplots(1, 2, figsize=(40, 40))
  3499. plot_order_fig1 = ['dwMRI', 'rsMRI']
  3501. for idx, group_name in enumerate(plot_order_fig1):
  3502. modality_list = modality_groups[group_name]
  3503. ax = axs1[idx]
  3505. # Filter data
  3506. group_data = df_individual[df_individual['Modality'].isin(modality_list)].copy()
  3508. # Sort by unique age variance
  3509. group_data = group_data.sort_values('Unique variance: brain/body', ascending=True)
  3511. # Extract data for plotting
  3512. y_labels = group_data['Modality'].values
  3513. y_labels = [label.replace('Full Correlation', 'FC').replace('Partial Correlation', 'PC').replace('Probabilistic', 'Prob.') for label in y_labels]
  3515. # Apply dictionary
  3516. rename_dict = {'Bone Densitometry (of heel)': 'Bone densitometry (heel)',
  3517. 'Abdominal Organ Composition by MRI': 'Abdominal organ composition (MRI)',
  3518. 'Abdominal Composition by MRI': 'Abdominal composition (MRI)',
  3519. 'Body Composition by DXA': 'Body composition (DXA)', 'Abdominal composition by MRI':
  3520. 'Abdominal Composition (MRI)',
  3521. 'Bone Size, Mineral and Density by DXA': 'BMD (DXA)',
  3522. 'Body Composition by Impedance': 'Body composition (impedance)' }
  3523. y_labels = [rename_dict.get(lbl, lbl) for lbl in y_labels]
  3525. y1 = group_data['Unique variance: age'].values*100
  3526. y2 = group_data['Common variance'].values*100
  3527. y3 = group_data['Unique variance: brain/body'].values*100
  3529. # Calculate percentages for annotations
  3530. total = y1 + y2 + y3
  3531. perc_u_age = (y1 / total * 100)
  3532. perc_common = (y2 / total * 100)
  3533. perc_u_brain_body = (y3 / total * 100)
  3534. marginal_u_age = group_data['Unique variance: age'].values
  3535. marginal_common = group_data['Common variance'].values
  3536. marginal_u_brain_body = group_data['Unique variance: brain/body'].values
  3538. # Get the proportion of g-age explained by brain/body
  3539. prop_age_exp_by_brain_body = group_data['g-age explained by brain/body'].values
  3541. bar_height = 1
  3543. # Create stacked horizontal bars
  3544. ax.barh(y_labels, y1, height=bar_height, color=color_unique_age, alpha=0.9, edgecolor='black', linewidth=0.9)
  3545. ax.barh(y_labels, y2, height=bar_height, left=y1, color=color_common, alpha=0.9, edgecolor='black', linewidth=0.9)
  3546. ax.barh(y_labels, y3, height=bar_height, left=y1+y2, color=modality_colors[group_name], alpha=0.9, edgecolor='black', linewidth=0.9)
  3548. # Add percentage annotations on each segment with font size dictionary
  3549. bar_annot_font_sizes = {
  3550. 'dwMRI': 30,
  3551. 'rsMRI': 40,
  3552. 'sMRI': 30,
  3553. 'Body': 30}
  3555. # Adjust thresholds
  3556. for i, (u_age, common, u_brain_body) in enumerate(zip(y1, y2, y3)):
  3557. if u_age > 0.001:
  3558. ax.text(u_age/2, i, f'{u_age:.1f}', ha='center', va='center',
  3559. color='black', fontsize=bar_annot_font_sizes[group_name])
  3560. if common > 0.001:
  3561. ax.text(u_age + common / 2, i, f'{common:.1f}', ha='center', va='center',
  3562. color='black', fontsize=bar_annot_font_sizes[group_name])
  3563. if u_brain_body:
  3564. pass #> 0.001:
  3565. #ax.text(u_age + common + u_brain_body/2, i, f'{m_u_brain_body:.2f}', ha='center', va='center', color='black', fontsize=bar_annot_font_sizes[group_name])
  3567. # Annotations
  3568. bbox_font_sizes = {
  3569. 'dwMRI': 38,
  3570. 'rsMRI': 48,
  3571. 'sMRI': 40,
  3572. 'Body': 40}
  3574. for i, prop_value in enumerate(prop_age_exp_by_brain_body):
  3575. prop_text = f'{prop_value:.1f}%'
  3576. # Position the text to the right of the plot
  3577. # transform=ax.get_yaxis_transform() to position relative to y-axis
  3578. ax.text(1.02, i, prop_text, ha='left', va='center',
  3579. fontsize=bbox_font_sizes[group_name], transform=ax.get_yaxis_transform(),
  3580. color='black', bbox=dict(boxstyle="round,pad=0.1", facecolor='white', alpha=0.7))
  3582. # Format the subplot
  3583. title_color = modality_colors[group_name]
  3584. ax.set_title(group_name,
  3585. fontsize=80, y=1,
  3586. bbox=dict(facecolor='white', alpha=0.9,
  3587. edgecolor=title_color, linewidth=5,
  3588. boxstyle='round,pad=0.3'))
  3589. y_labelsize_dict = {
  3590. 'dwMRI': 38,
  3591. 'rsMRI': 60,
  3592. 'sMRI': 60,
  3593. 'Body': 60}
  3595. ax.tick_params(axis='y', labelsize=y_labelsize_dict[group_name], length=0)
  3596. ax.tick_params(axis='x', labelsize=65)
  3597. ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%.1f'))
  3599. # Set x-axis limit
  3600. max_total = max(y1 + y2 + y3)
  3601. ax.set_xlim(0, max_total)
  3602. val1 = max_total / 3
  3603. val2 = 2 * max_total / 3
  3604. ax.xaxis.set_major_locator(mticker.FixedLocator([0, val1, val2, max_total]))
  3606. # Remove spines
  3607. for spine in ['top', 'right', 'left']:
  3608. ax.spines[spine].set_visible(False)
  3610. # Add common x axis label
  3611. fig1.text(0.57, 0.08, 'Variance of Cognition Explained (%)',
  3612. ha='center', fontsize=80)
  3614. # Add common title
  3615. fig1.suptitle('Commonality Analysis\nbetween Age and each Body and Brain Phenotype\nin Explaining Cognition (%)',
  3616. fontsize=100, y=0.98, x=0.55)
  3618. # Adjust layout
  3619. plt.tight_layout(rect=[0, 0.12, 1, 0.95])
  3620. plt.subplots_adjust(wspace=1.6, hspace=0.2)
  3622. # transform
  3623. save_fig(['png', 'svg'], os.path.join(fig_path, 'Fig7_part1_dwMRI_rsMRI_barplots'))
  3624. plt.show()
  3626. # %% [markdown]
  3627. # ## Figure 7b
  3629. # %%
  3630. # Bottom Row - sMRI and Body
  3631. fig2, axs2 = plt.subplots(1, 2, figsize=(55, 30))
  3632. axs2 = axs2.flatten()
  3633. plot_order_fig2 = ['sMRI', 'Body']
  3635. for idx, group_name in enumerate(plot_order_fig2):
  3636. modality_list = modality_groups[group_name]
  3637. ax = axs2[idx]
  3639. # Filter data
  3640. group_data = df_individual[df_individual['Modality'].isin(modality_list)].copy()
  3642. # Sort by unique age variance
  3643. group_data = group_data.sort_values('Unique variance: brain/body', ascending=True)
  3645. # Extract data for plotting
  3646. y_labels = group_data['Modality'].values
  3647. y_labels = [label.replace('Full Correlation', 'FC').replace('Partial Correlation', 'PC').replace('Probabilistic', 'Prob.') for label in y_labels]
  3649. # Apply dictionary
  3650. rename_dict = {'Bone Densitometry (of heel)': 'Bone densitometry (heel)',
  3651. 'Abdominal Organ Composition by MRI': 'Abdominal organ composition (MRI)',
  3652. 'Abdominal Composition by MRI': 'Abdominal composition (MRI)',
  3653. 'Body Composition by DXA': 'Body composition (DXA)', 'Abdominal composition by MRI':
  3654. 'Abdominal Composition (MRI)',
  3655. 'Bone Size, Mineral and Density by DXA': 'BMD (DXA)',
  3656. 'Body Composition by Impedance': 'Body composition (impedance)',
  3657. 'Desikan Grey/white Matter Contrast': 'Desikan Grey/White Matter Contrast',
  3658. 'Pulse Wave Analysis' : 'Pulse wave analysis',
  3659. 'Arterial Stiffness' : 'Arterial stiffness',
  3660. 'Renal & Hepatic' : 'Renal & hepatic',
  3661. 'ECG at Rest' : 'ECG at rest',
  3662. 'Bone Densitometry of heel' : 'Bone densitometry of heel',
  3663. 'Carotid Ultrasound' : 'Carotid ultrasound'}
  3665. y_labels = [rename_dict.get(lbl, lbl) for lbl in y_labels]
  3667. y1 = group_data['Unique variance: age'].values*100
  3668. y2 = group_data['Common variance'].values*100
  3669. y3 = group_data['Unique variance: brain/body'].values*100
  3671. # Calculate percentages for annotations
  3672. total = y1 + y2 + y3
  3673. perc_u_age = (y1 / total * 100)
  3674. perc_common = (y2 / total * 100)
  3675. perc_u_brain_body = (y3 / total * 100)
  3677. # Get the proportion of g-age explained by brain/body
  3678. prop_age_exp_by_brain_body = group_data['g-age explained by brain/body'].values
  3680. bar_height = 1
  3682. # Create stacked horizontal bars
  3683. ax.barh(y_labels, y1, height=bar_height, color=color_unique_age, alpha=0.9, edgecolor='black', linewidth=0.9)
  3684. ax.barh(y_labels, y2, height=bar_height, left=y1, color=color_common, alpha=0.9, edgecolor='black', linewidth=0.9)
  3685. ax.barh(y_labels, y3, height=bar_height, left=y1+y2, color=modality_colors[group_name], alpha=0.9, edgecolor='black', linewidth=0.9)
  3687. # Add percentage annotations on each segment with font size dictionary
  3688. bar_annot_font_sizes = {
  3689. 'dwMRI': 30,
  3690. 'rsMRI': 40,
  3691. 'sMRI': 40,
  3692. 'Body': 40}
  3694. # Adjust thresholds
  3695. for i, (u_age, common, u_brain_body) in enumerate(zip(y1, y2, y3)):
  3696. if u_age > 0.001:
  3697. ax.text(u_age/2, i, f'{u_age:.1f}', ha='center', va='center',
  3698. color='black', fontsize=bar_annot_font_sizes[group_name])
  3699. if common > 0.001:
  3700. ax.text(u_age + common/2, i, f'{common:.1f}', ha='center', va='center',
  3701. color='black', fontsize=bar_annot_font_sizes[group_name])
  3702. if u_brain_body:
  3703. pass #> 0.001:
  3704. #ax.text(u_age + common + u_brain_body/2, i, f'{m_u_brain_body:.2f}', ha='center', va='center', color='black', fontsize=bar_annot_font_sizes[group_name])
  3706. # Annotations
  3707. bbox_font_sizes = {
  3708. 'dwMRI': 38,
  3709. 'rsMRI': 48,
  3710. 'sMRI': 48,
  3711. 'Body': 48}
  3714. for i, prop_value in enumerate(prop_age_exp_by_brain_body):
  3715. prop_text = f'{prop_value:.1f}%'
  3716. # Position the text to the right of the plot
  3717. # transform=ax.get_yaxis_transform() to position relative to y-axis
  3718. ax.text(1.02, i, prop_text, ha='left', va='center',
  3719. fontsize=bbox_font_sizes[group_name], transform=ax.get_yaxis_transform(),
  3720. color='black', bbox=dict(boxstyle="round,pad=0.1", facecolor='white', alpha=0.7))
  3722. # Format the subplot
  3723. title_color = modality_colors[group_name] # pick edge color from dict
  3724. ax.set_title(group_name,
  3725. fontsize=80, y=1,
  3726. bbox=dict(facecolor='white', alpha=0.9,
  3727. edgecolor=title_color, linewidth=5,
  3728. boxstyle='round,pad=0.3'))
  3729. y_labelsize_dict = {
  3730. 'dwMRI': 38,
  3731. 'rsMRI': 60,
  3732. 'sMRI': 60,
  3733. 'Body': 60}
  3735. ax.tick_params(axis='y', labelsize=y_labelsize_dict[group_name], length=0)
  3736. ax.tick_params(axis='x', labelsize=65)
  3737. ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%.1f'))
  3739. # Set x-axis limit - extend slightly to accommodate the right-side annotations
  3740. max_total = max(y1 + y2 + y3)
  3741. ax.set_xlim(0, max_total)
  3742. val1 = max_total / 3
  3743. val2 = 2 * max_total / 3
  3744. ax.xaxis.set_major_locator(mticker.FixedLocator([0, val1, val2, max_total]))
  3746. # Remove spines
  3747. for spine in ['top', 'right', 'left']:
  3748. ax.spines[spine].set_visible(False)
  3750. # Add common x axis label
  3751. fig2.text(0.63, 0.08, 'Variance of Cognition Explained (%)',
  3752. ha='center', fontsize=80)
  3754. # Add circle legend
  3755. legend_elements_circles = [
  3756. plt.Circle((0,0), 0.5, facecolor=color_unique_age, alpha=0.9, label='Unique Variance: Age'),
  3757. plt.Circle((0,0), 0.5, facecolor=color_common, alpha=0.9, label='Common Variance'),
  3758. plt.Circle((0,0), 0.5, facecolor=modality_colors['sMRI'], alpha=0.9, label='Unique Variance: sMRI'),
  3759. plt.Circle((0,0), 0.5, facecolor=modality_colors['dwMRI'], alpha=0.9, label='Unique Variance: dwMRI'),
  3760. plt.Circle((0,0), 0.5, facecolor=modality_colors['rsMRI'], alpha=0.9, label='Unique Variance: rsMRI'),
  3761. plt.Circle((0,0), 0.5, facecolor=modality_colors['Body'], alpha=0.9, label='Unique Variance: Body'),
  3762. ]
  3764. leg2_top = fig2.legend(handles=legend_elements_circles,
  3765. loc='lower center', ncol=3,
  3766. fontsize=75, frameon=False, fancybox=False,
  3767. bbox_to_anchor=(0.6, -0.15), # Higher
  3768. handler_map={plt.Circle: HandlerPatch(patch_func=make_circle_legend)},
  3769. handletextpad=1.0, handleheight=1.0,
  3770. columnspacing=1.2, labelspacing=1.2, borderpad=0.5)
  3772. # Add rectangle legend
  3773. legend_element_rectangle = [
  3774. plt.Rectangle((0,0), 1, 1, facecolor='none', edgecolor='black', lw=2,
  3775. label='Age-related cognitive variance explained by each body and brain phenotype (%)')
  3776. ]
  3778. leg2_bottom = fig2.legend(handles=legend_element_rectangle,
  3779. loc='lower center', ncol=1,
  3780. fontsize=75, frameon=False, fancybox=False,
  3781. bbox_to_anchor=(0.55, -0.25), # Lower than circles
  3782. handletextpad=0.5, handleheight=1.0,
  3783. columnspacing=1.2, labelspacing=1.2, borderpad=0.5)
  3785. # Adjust layout
  3786. plt.tight_layout(rect=[0, 0.12, 1, 0.95])
  3787. plt.subplots_adjust(wspace=3, hspace=0.2)
  3789. # Save
  3790. save_fig(['png', 'svg'], os.path.join(fig_path, 'Fig7_part2_sMRI_Body_barplots'))
  3791. plt.show()

1Figures.ipynb at commit 54c5c9d, under MIT · at the source

Overview

Authors: Irina Buianova1, Narun Pat1
  1. Department of Psychology, University of Otago, Dunedin, New Zealand
Institutions: University of Otago (New Zealand)
Journal: npj aging, volume 12, issue 1, article 108
Dates: received 18 March 2026; accepted 17 July 2026; published online 31 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41514-026-00456-9 · PMID 42581321 · PMCID PMC13462990 · OpenAlex W7171937603
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), cognitive (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Physiology & signal measures, fMRI & imaging, Preprocessing
Keywords: Neurology, Neuroscience
Topic: Dementia and Cognitive Impairment Research (Psychiatry and Mental health, Medicine), according to OpenAlex
Funding: Neurological Foundation of New Zealand (2350 PRG); The Ministry of Business, Innovation and Employment (UOA2421 and RTVU2403); Health Research Council of New Zealand (21/618 and 24/838)
Citations: not cited yet (Europe PMC); 197 references in the paper

Abstract

Epidemiological links between cognition and body physiology in aging are well established, but their strength and drivers remain unclear. Which physiological systems – from body composition to cardiovascular, pulmonary, renal, hepatic, immune, metabolic, and musculoskeletal – best predict cognition, and to what extent are cognition–body associations linked to brain variation across aging? We examined 19 physiological phenotypes alongside three neuroimaging modalities in over 30,000 UK Biobank participants. Machine learning models integrating body measures predicted cognition at r = 0.4, demonstrating a cognition–body covariation at 16%. Body composition and bone health emerged as the strongest predictors. Notably, 85.1% of cognition–body covariance overlapped with neuroimaging, especially white matter features. Moreover, 71.7% of cognition–age covariance was jointly shared with neuroimaging and physiology, and 96.8% was shared with either brain or body markers, or their overlap. Together, these findings clarify how body physiology and brain structure and function covary with cognitive aging.

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

HAM-lab-Otago-University/UKBiobank-Brain-Body

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 54c5c9ddf20dfe52285a9e50803b992724bd7c88, 2 June 2026
Languages: Jupyter (18), Python (18), R (3)
Size: 60 files, 39 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 21 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (36 files), pandas (36 files), scikit-learn (32 files), SciPy (25 files), XGBoost (12 files), Matplotlib (11 files), seaborn (10 files), Nilearn (5 files), statsmodels (4 files), lavaan (3 files), data.table (2 files), ggplot2 (2 files), psych (2 files), NiBabel (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
41 files
At the source: github.com/HAM-lab-Otago-University/UKBiobank-Brain-Body

Code availability

All analysis code is available at https://github.com/HAM-lab-Otago-University/UKBiobank-Brain-Body. All analyses can be reproduced using the provided code and UK Biobank data obtained through the appropriate access procedures.

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

Data availability

The authors are not permitted to publicly share individual-level UK Biobank data due to legal and ethical restrictions. Access requests should be directed to the UK Biobank (https://www.ukbiobank.ac.uk/), subject to its data access policies and approval procedures. This manuscript is a computational study and did not generate new primary data.

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, 2 authors, 2 keywords, 3 funders, 190 references.

Cite

This paper

Buianova, I., & Pat, N. (2026). Exploring the link between body physiology and cognition: the role of the brain and aging. npj aging, 12(1), 108. https://doi.org/10.1038/s41514-026-00456-9

BibTeX

@article{buianova2026exploring,
author = {Buianova, Irina and Pat, Narun},
title = {{Exploring the link between body physiology and cognition: the role of the brain and aging}},
journal = {npj aging},
year = {2026},
month = jul,
volume = {12},
number = {1},
pages = {108},
publisher = {Nature Publishing Group},
issn = {2731-6068},
doi = {10.1038/s41514-026-00456-9},
url = {https://doi.org/10.1038/s41514-026-00456-9},
pmid = {42581321},
pmcid = {PMC13462990}
}

RIS

TY - JOUR
AU - Buianova, Irina
AU - Pat, Narun
TI - Exploring the link between body physiology and cognition: the role of the brain and aging
T2 - npj aging
J2 - NPJ Aging
PY - 2026
DA - 2026/07/31
VL - 12
IS - 1
SP - 108
SN - 2731-6068
PB - Nature Publishing Group
DO - 10.1038/s41514-026-00456-9
UR - https://doi.org/10.1038/s41514-026-00456-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41514-026-00456-9",
"type": "article-journal",
"title": "Exploring the link between body physiology and cognition: the role of the brain and aging",
"container-title": "npj aging",
"author": [
{
"family": "Buianova",
"given": "Irina"
},
{
"family": "Pat",
"given": "Narun"
}
],
"container-title-short": "NPJ Aging",
"volume": "12",
"issue": "1",
"page": "108",
"DOI": "10.1038/s41514-026-00456-9",
"PMID": "42581321",
"PMCID": "PMC13462990",
"ISSN": "2731-6068",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41514-026-00456-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
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.7554/elife.108109 [code]
Multimodal MRI marker of cognition explains the association between cognition and mental health in the UK Biobank.
Journal: eLife
In common: lavaan, psych, data.table, 7 other tools, cognitive, 34 references
[2] doi:10.1038/s41467-026-73668-y [code]
Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: Nilearn, NiBabel, ggplot2, 7 other tools, 10 references
[3] doi:10.1038/s41467-026-75661-x [code]
A neural signature of sleep deprivation in the human brain.
Journal: Nature communications
In common: Nilearn, NiBabel, statsmodels, 6 other tools, cognitive, 7 references
[4] doi:10.1162/imag.a.1282 [code]
Metabolic syndrome severity and the energetic cost of brain network transitions: A normative modeling study of accelerated brain aging.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, NiBabel, statsmodels, 6 other tools, 6 references
[5] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: XGBoost, psych, NiBabel, 9 other tools, 3 references
[6] doi:10.1162/imag.a.105 [code]
Right posterior theta reflects human parahippocampal phase resetting by salient cues during goal-directed navigation
Journal: n/a
In common: lavaan, psych, Nilearn, 9 other tools, cognitive, 1 reference
[7] doi:10.1038/s42003-026-09956-6 [code]
Linking changes in sulcal morphometry to cognitive development from childhood to adolescence.
Journal: Communications biology
In common: psych, Nilearn, NiBabel, 9 other tools, 2 references
[8] doi:10.21203/rs.3.rs-9914920/v1 [code]
Prediction of cognitive performance by demographics, sleep, and brain morphometry: machine learning findings from ENIGMA-Sleep Working Group
Journal: Research Square (preprint)
In common: XGBoost, Nilearn, NiBabel, 7 other tools, 2 references
[9] 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: lavaan, psych, statsmodels, 7 other tools, 2 references
[10] doi:10.7554/elife.103097 [code]
Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: Nilearn, data.table, NiBabel, 7 other tools, 3 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.