Metabolism-weighted brain connectome reveals synaptic integration and vulnerability to neurodegeneration.

Code ↔ Paper

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

The 12 matches

[1] § Materials and Methods › Weighted Degree Centrality (wDC). ↔ src/Functions.ipynb, lines 502–595 · score 0.66 · Pearson correlation, edge weighted, FC matrix, thresholded, masked, SC
[2] § Results › Metabolism-Weighted Hubness Reflects Synaptic Gene Expression and Susceptibility to Neurodegenerative Diseases. ↔ scripts/Figures_codes.ipynb, lines 2218–2277 · score 0.62 · oxidative phosphorylation, Alzheimer, neurodegenerative, KEGG, diseases, pathways
[3] § Materials and Methods › Participants. ↔ scripts/Figures_codes.ipynb, lines 97–209 · score 0.61 · Vienna.rep, TUM.rep, split
[4] § Materials and Methods › Control Analyses. ↔ scripts/Figures_codes.ipynb, lines 3759–3889 · score 0.60 · MI threshold, SC threshold, sensitivity, voxels, energy, PET
[5] § Results › MwC Reveals Energy-Based Architecture of Cortical Integration. ↔ scripts/Figures_codes.ipynb, lines 3891–4033 · score 0.59 · Julich Brain atlas, spatial correlation, CMRglc, MNI, S6, parcellation
[6] § Materials and Methods › Statistical Analyses. ↔ src/Functions.ipynb, lines 1530–1558 · score 0.57 · correlation coefficients, dependent correlations, Steiger, variable, scored
[7] § Materials and Methods › Statistical Analyses. ↔ src/functions.py, lines 732–759 · score 0.57 · correlation coefficients, dependent correlations, Steiger, variable, scored
[8] § Results › Metabolism-Weighted Hubs Align with Cognitive Network Hierarchies. ↔ scripts/Figures_codes.ipynb, lines 1474–1489 · score 0.57 · eye movements, visuospatial, inhibition, Neurosynth, memory, cognitive
[9] § Materials and Methods › Participants. ↔ scripts/MwC_generation.ipynb, lines 18–61 · score 0.56 · Vienna.rep, TUM.rep
[10] § Materials and Methods › Model Description. ↔ src/Functions.ipynb, lines 912–915 · score 0.54 · Gaussian Copula Mutual, mutual information
[11] § Materials and Methods › Model Description. ↔ src/functions.py, lines 336–338 · score 0.54 · Gaussian Copula Mutual, mutual information
[12] § Materials and Methods › Weighted Degree Centrality (wDC). ↔ src/functions.py, lines 922–984 · score 0.51 · edge weighted, FC matrix, thresholded, masked, SC, graph

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 paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 4,803 lines · 159 KB · MIT · 5 matches

# %% [markdown]
# # Setting
# %% [markdown]
# ## libraries
# %%
%load_ext autoreload
%autoreload 2
import sys
from pathlib import Path
SRC = (Path.cwd().parent / "src").resolve()
sys.path.insert(0, str(SRC))
import importlib
import functions
importlib.reload(functions)
from functions import *
# %%
import os
import os,sys
import numpy as np
import scipy as sp
import igraph as ig
import scipy.speciali ws
import scipy.stats as stats
from scipy.stats import norm
np.float = float
np.int = int
import pandas as pd
import pickle
import csv
import shutil
import time
import warnings
import nilearn as ni
from nilearn import datasets, plotting #, surface, plotting, input_data
from nilearn.image import math_img, threshold_img, smooth_img
from nilearn.plotting import cm as nip_cm
import nibabel as nib
import networkx as nx
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, Normalize
from matplotlib.patches import Patch
from matplotlib.cm import ScalarMappable
from scipy.stats import pearsonr, zscore, spearmanr
from scipy.spatial.distance import euclidean
from scipy.spatial import distance
from scipy import linalg, optimize
import matplotlib.colorbar as colorbar
from enigmatoolbox.permutation_testing import rotate_parcellation, perm_sphere_p
import matplotlib as mpl
# import plotly.io as pio
np.seterr(invalid='ignore')
#import enigmatoolbox
from enigmatoolbox.utils.parcellation import surface_to_parcel, parcel_to_surface
import statsmodels.api as sm
import plotly.graph_objects as go
# from brainsmash.mapgen.base import Base
# from brainsmash.mapgen.eval import base_fit
# from brainsmash.mapgen.stats import nonparp, pairwise_r
# import plotly.graph_objects as go
# from neuromaps.datasets import fetch_annotation
# import matplotlib.ticker as ticker
# import textwrap
# # Define FSL directories
os.environ["FSLDIR"] = "/usr/share/fsl/5.0"
os.environ["FSLOUTPUTTYPE"] = "NIFTI_GZ"
os.environ["FSLTCLSH"] = "/usr/bin/tclsh"
os.environ["FSLWISH"] = "/usr/bin/wish"
os.environ["FSLMULTIFILEQUIT"] = "True"
os.environ["LD_LIBRARY_PATH"] = "/usr/share/fsl/5.0:/usr/lib/fsl/5.0"
# Fix deprecated NumPy aliases
np.float = float
np.int = int
# %% [markdown]
# ## IMPORTANT: select the data (main or external) AND including/excluding the limbic nework
# %%
# ================================================
# 1. session selection (main or external dataset):
# ================================================
session = "AUF" # AUF(main), ZU(TUM.rep), m1 (Vienna.rep session 1) or m2 (Vienna.rep session 2)
#DIR = f'/data/raw/{session}/'
nrois = 360
# ================================================
# 2. select if you want to include limbic or not:
# ================================================
# including limbic:
# limb_ind = []
# LIMB = 'with'
# excluding limbic:
limb_ind = np.array([88, 90, 92, 93, 110, 118, 120, 122, 131, 135, 165, 166, 172, 268, 270, 272, 273, 290, 298, 300, 302, 311, 314, 315, 346, 352])-1
LIMB = 'without'
# ================================================
# 3. Parameter assignments:
# ================================================
# sub = os.listdir(DIR)
# sub = sorted([int(s.split('-')[1]) for s in sub if s.startswith('sub-') and s.split('-')[1].isdigit()])
# if session == "ZU":
# sub.remove(14)
# if session == 'AUF':
# sub = [3,7,12,14,17,20,23,25,26,28,29,30,31,32,33,35,36,37,38]
# sub_size = len(sub)
nrois_rem = nrois - len(limb_ind)
net_label = pd.read_csv(f'../data/external/mmp2yeo7nw_mapping.csv')
net_label = net_label.drop(limb_ind,axis=0)
net_names = pd.unique(net_label['yeo_7_nw'])
net_label.reset_index(drop=True, inplace=True)
network_to_number = {network: i + 1 for i, network in enumerate(net_label['yeo_7_nw'].unique())}
net_label['network_number'] = net_label['yeo_7_nw'].map(network_to_number)
## network colors assignment:
color = [(139, 19, 140),(1,131, 182),(51, 116, 32),
(222, 75, 82),(226, 55, 255),(239, 156, 60), (255, 254, 211)]
net_color = [(r/255, g/255, b/255 , 1) for r,g,b in color]
customPalette = sns.set_palette(sns.color_palette(net_color))
xtik = ['Vis','Som','Dors','Def', 'Sal' , 'Cont', 'Limbic']
net_names = ['Vis', 'SomMot' ,'DorsAttn' ,'Default' ,'SalVentAttn', 'Cont', 'Limbic' ]
net_num = net_label["network_number"]
new_net_num = np.tile(net_num.transpose(), (nrois_rem, 1))
num_net = len(net_names)
# ================================================
# 4. Parcellation setting:
# ================================================
mmp = nib.load(f'../data/external/MMP_in_MNI_corr_3mm.nii.gz').get_fdata()
v11 , v22 , v33 = mmp.shape
nvox2 = v11 * v22 * v33
mmp_re = mmp.reshape((nvox2, 1),order='F')
regions = np.unique(mmp_re)[1:]
regions_rem = regions.copy()
regions_rem = np.delete(regions_rem , limb_ind)
# ================================================
# 4. figure's properties:
# ================================================
FONT = 18
COLOR = "#8DD084" #"#97A8C0" #"#8DD084"
IC_PET_COLf = '#70DEB6'
IC_PET_COLe = '#1FC488'
organg = [255/255,164/255,0/255]
GRAY = "#7F7F7F"
deg_color = "#97A8C0"
IC_color = "#8DD084"
approximate_colors = [
"#341C54",
"#3E2166",
"#482677", # Dark Purple
"#44357F",
"#404387", # Purple-Blue
"#33638D", # Blue
"#2A788E", # Blue-Green
"#1F9E89", # Greenish-Blue
"#35B779", # Green
"#6DCD59", # Lime Green
"#B4DE2C", # Yellow-Green
"#D9E329",
"#EBE527",
"#FDE725" # Yellow
]
approximate_cmap = LinearSegmentedColormap.from_list("approximate_colormap", approximate_colors, N=256)
CMAP_IC = approximate_cmap
cool_cmap = plt.get_cmap('cool')
#cool_cmap = cm.get_cmap('cool')
cool_colors = cool_cmap(np.linspace(0.15, 0.95, 256))
CMAP_DC = LinearSegmentedColormap.from_list("adjusted_cool", cool_colors)
# %% [markdown]
# ## color maps/ colors
# %%
approximate_colors = [
"#341C54",
"#3E2166",
"#482677", # Dark Purple
"#44357F",
"#404387", # Purple-Blue
"#33638D", # Blue
"#2A788E", # Blue-Green
"#1F9E89", # Greenish-Blue
"#35B779", # Green
"#6DCD59", # Lime Green
"#B4DE2C", # Yellow-Green
"#D9E329",
"#EBE527",
"#FDE725" # Yellow
]
approximate_cmap = LinearSegmentedColormap.from_list("approximate_colormap", approximate_colors, N=256)
CMAP_IC = approximate_cmap
import matplotlib.pyplot as plt
# Now this works:
cool_cmap = plt.get_cmap('cool')
#cool_cmap = cm.get_cmap('cool')
cool_colors = cool_cmap(np.linspace(0.15, 0.95, 256))
CMAP_DC = LinearSegmentedColormap.from_list("adjusted_cool", cool_colors)
network_colors = [(139, 19, 140), (1, 131, 182), (51, 116, 32), (222, 75, 82), (226, 55, 255), (239, 156, 60), (255, 254, 211)]
network_colors = [(r / 255, g / 255, b / 255) for (r, g, b) in network_colors]
network_colors = ['#9A199A', '#45B7F0','#58B73B','#DE4B52','#E237FF','#F29A36','#FEFFD3']
lighter_colors = [sns.light_palette(color, n_colors=100, input="hex")[50] for color in network_colors]
#lighter_colors2 = [(*sns.color_palette([color])[0], 0.01) for color in lighter_colors] # Adding 50% opacity
lighter_colors2 = ['#C788C7' , '#C0EBFF' , '#BEE6B2' , '#F2BABD', '#F6C2FF','#F7C58D']
network_order_original = ['Vis', 'SomMot', 'DorsAttn', 'Default', 'SalVentAttn', 'Cont']
network_color_map = dict(zip(network_order_original, network_colors))
lighter_color_map = dict(zip(network_order_original, lighter_colors))
import matplotlib.colors as mcolors
def matplotlib_to_plotly(cmap, n=256):
"""Convert a Matplotlib colormap to a Plotly-friendly colorscale."""
return [
[i / (n - 1), mcolors.rgb2hex(cmap(i / (n - 1))[:3])]
for i in range(n)
]
# %% [markdown]
# # Preparing data
# %% [markdown]
# ## Loading MwC and DC matrices
# %%
'''
The following files were generated in MwC_generation.ipynb notebook:
1. matrices_data_{session}_{LIMB}_limbic.pkl
2. DATA_avg_wSC_{session}_{LIMB}_limb.csv
3. data_single_sub_wSC_{session}_{LIMB}_limb.pkl
4. DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv
5. data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl
'''
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
with open(file_path, 'rb') as file:
data_input = pickle.load(file)
num_voxels = data_input['num_voxels']
MIconn_allsub = data_input['MIconn_allsub']
SCconn_allsub = data_input['SCconn_allsub']
Pearconn_allsub = data_input['Pearconn_allsub']
edata_medianallsub = data_input['edata_medianallsub']
#--------------------------------------------------------------------------------------------
#-------------- Reading MWC matrices --------------------------
#--------------------------------------------------------------------------------------------
## with structural connectivity masking:
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_wSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl'
## without structural connectivity masking:
# data_avg_pr_path = f'../data/processed/DATA_avg_woSC_pr_{session}_{LIMB}_limb.csv'
# data_sub_pr_path = f'../data/processed/data_single_sub_woSC_pr_{session}_{LIMB}_limb.pkl'
# data_avg_path = f'../data/processed/DATA_avg_woSC_{session}_{LIMB}_limb.csv'
# data_sub_path = f'../data/processed/data_single_sub_woSC_{session}_{LIMB}_limb.pkl'
DATA_avg = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub = pickle.load(file)
ICallsub_w = data_single_sub['ICallsub_w']
ICallsub_b = data_single_sub['ICallsub_b']
degallsub_w = data_single_sub['degallsub_w']
degallsub_b = data_single_sub['degallsub_b']
AvgMIallsub = data_single_sub['AvgMIallsub']
#--------------------------------------------------------------------------------------------
#-------------- Reading DC matrices ----------------
#--------------------------------------------------------------------------------------------
DATA_avg_pr = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr = pickle.load(file)
ICallsub_w_pr = data_single_sub_pr['ICallsub_w']
ICallsub_b_pr = data_single_sub_pr['ICallsub_b']
degallsub_w_pr = data_single_sub_pr['degallsub_w']
degallsub_b_pr = data_single_sub_pr['degallsub_b']
AvgMIallsub_pr = data_single_sub_pr['AvgMIallsub']
sub_size = ICallsub_w_pr.shape[1]
# %% [markdown]
# ## Correlation values
# %%
from scipy.stats import pearsonr
corr_IC_E_allsub = np.zeros((sub_size))
corr_btw_E_allsub = np.zeros((sub_size))
corr_eig_E_allsub = np.zeros((sub_size))
corr_deg_w_E_allsub = np.zeros((sub_size))
corr_IC_E_net_allsub = np.zeros((sub_size, len(net_names)))
p_IC_E_allsub = np.zeros((sub_size))
p_deg_w_E_allsub = np.zeros((sub_size))
p_btw_E_allsub = np.zeros((sub_size, 1))
p_eig_E_allsub = np.zeros((sub_size, 1))
p_IC_E_net_allsub = np.zeros((sub_size, len(net_names)))
for j in range(sub_size):
degallsub_w_z = zscore(degallsub_w[:, j] , nan_policy='omit')
corr_IC_E_allsub[j], p_IC_E_allsub[j] = pcor(ICallsub_w[:, j], edata_medianallsub[:, j])
corr_deg_w_E_allsub[j], p_deg_w_E_allsub[j] = pcor(degallsub_w_z, edata_medianallsub[:, j])
corr_IC_w_E, p_IC_w_E = pcor(DATA_avg.ICallsub_w_avg, DATA_avg.pet_avg)
corr_IC_b_E, p_IC_b_E = pcor(DATA_avg.ICallsub_b_avg, DATA_avg.pet_avg)
corr_IC_w_deg, p_IC_w_deg = pcor(DATA_avg.ICallsub_w_avg, DATA_avg.degallsub_w_avg)
corr_IC_b_deg, p_IC_b_deg = pcor(DATA_avg.ICallsub_b_avg, DATA_avg.degallsub_b_avg)
corr_deg_w_E, p_deg_w_E = pcor(zscore(DATA_avg.degallsub_w_avg , nan_policy='omit'), DATA_avg.pet_avg)
# %% [markdown]
# # Figures
# %% [markdown]
# ## Fig 2a: MwC vs CMRglc
# %%
#--------------------------------------------------------------------------------------------
#---------------------- spatial autocorrelation (SAC) -----------------------------------
#--------------------------------------------------------------------------------------------
##### BRAIN SMASH:
# niter = 10
# test_stat_IC ,surrogate_brainmap_corrs_IC, sa_corrected_p_value_IC, spatially_naive_p_value_IC = Spatial_AC(DATA_avg , "IC", LIMB,niter)
# sac1 = '#0199DD'
# sac2 = '#8CCF83'
# sac3 = '#FEAD01'
# fig, ax = plt.subplots(figsize=(3, 7))
# plt.grid(True)
# g = sns.kdeplot(surrogate_brainmap_corrs_IC, color= "#E0E0E0", ax=ax, fill= True, linewidth=2.5 )
# ax.axvline(corr_IC_w_E, 0, 0.9, color= COLOR, linestyle='dashed', lw=3)
# ax.set_xticks(np.arange(-1, 1.1, 0.5))
# ax.set_ylim(0, 2)
# ax.set_xlim(-0.8, 0.8)
# spine_width = 1.6
# spine_color = 'gray'
# for spine in ax.spines.values():
# spine.set_linewidth(spine_width)
# spine.set_color(spine_color)
# plt.tick_params(axis='both', which='major', labelsize= FONT)
# plt.xlabel('Correlation', fontsize = FONT)
# plt.ylabel('Density', fontsize = FONT)
# plt.savefig(f'../results/Figures/CMRglc_IC_hist.png', dpi=300, bbox_inches='tight')
# plt.show()
# lower_bound = np.percentile(surrogate_brainmap_corrs_IC, 5)
# upper_bound = np.percentile(surrogate_brainmap_corrs_IC, 95)
# print(f"95% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
### SPIN TEST###################
N_ROT = 10000
RNG_SEED = 42
session = 'AUF'
LIMB = 'with'
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
DATA_avg_spin = pd.read_csv(data_avg_path)
cmrglc = DATA_avg_spin['pet_avg'].to_numpy(dtype=float).ravel()
ic = DATA_avg_spin['ICallsub_w_avg'].to_numpy(dtype=float).ravel()
#rho_ic, p_spin_ic, null_ic = run_spin_test(ic, cmrglc, N_ROT,RNG_SEED)
FONT = 22
COLOR = IC_color
null_r = np.asarray(null_ic, dtype=float)
null_r = null_r[np.isfinite(null_r)]
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
sns.kdeplot(
null_r,
color="#E0E0E0",
ax=ax,
fill=True,
linewidth=2.5
)
ax.axvline(rho_ic, 0, 0.9, color=COLOR, linestyle='dashed', lw=3)
ax.set_xticks(np.arange(-1, 1.1, 0.5))
ax.set_ylim(0, 4)
ax.set_xlim(-0.8, 0.8)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.tick_params(axis='both', which='major', labelsize=FONT)
plt.xlabel('Correlation', fontsize=FONT)
plt.ylabel('Density', fontsize=FONT)
#plt.savefig('../results/Figures/CMRglc_IC_hist.png', dpi=300, bbox_inches='tight')
plt.show()
lower_bound = np.percentile(null_r, 5)
upper_bound = np.percentile(null_r, 95)
print(f"90% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
#--------------------------------------------------------------------------------------------
#------------------ scatter IC_avg , CMRglc_avg / group analysis ----------------------------
#--------------------------------------------------------------------------------------------
COLOR = IC_color
plt.figure(figsize=(6,4))
g = sns.jointplot(data=DATA_avg, x='pet_avg', y='ICallsub_w_avg', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
p_IC_w_E_smash = p_IC_w_E
if p_IC_w_E_smash < 0.0001:
title = "r = {:.2f} , p-smash < {:.3f}".format(float(corr_IC_w_E), float(0.001))
else:
title = "r = {:.2f} , p-smash = {:.3f}".format(float(corr_IC_w_E), float(p_IC_w_E_smash))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'MwC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Target Energy (CMRglc)', fontsize = FONT)
plt.ylabel('MwC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
#plt.savefig(f'../results/Figures/CMRglc_IC_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------ Box plot individual correlation between IC and E -----------------------
#--------------------------------------------------------------------------------------------
data = pd.DataFrame({'r': corr_IC_E_allsub, 'p': p_IC_E_allsub})
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
violin = sns.violinplot(data=data, y='r', ax=ax, color = COLOR)
violin.collections[0].set_facecolor('none')
sns.stripplot(data=data[data['p'] < 0.05], y='r', color=COLOR , size = 8, ax=ax)
sns.stripplot(data=data[data['p'] >= 0.05], y='r', color='red' , size = 8, ax=ax)
ax.set_ylabel('Correlation MwC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.15,0.72])
plt.tick_params(axis='both', which='major', labelsize= FONT)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
#plt.savefig(f'../results/Figures/CMRglc_IC_box.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- IC map on the brain surface -----------------------------------
#--------------------------------------------------------------------------------------------
IC = np.array(DATA_avg.ICallsub_w_avg).reshape(-1, 1)
IC_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[limb_ind] = False
ind = ind[mask]
IC_allrois[ind] = np.abs(IC)
fig = plt.figure(figsize=(5, 5),dpi=300)
fsaverage = datasets.fetch_surf_fsaverage()
V_min = IC.min()
V_max = IC.max()
V_max = 34
col_fsa = parcel_to_surface(IC_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP_IC,colorbar=False)
#plt.savefig(f'../results/Figures/IC_surf_l_lateral.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(IC_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data = False,
darkness=0.5,cmap=CMAP_IC,colorbar=False)
#plt.savefig(f'../results/Figures/IC_surf_l_medial.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(IC_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map = col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_right'], bg_on_data=False,
darkness=0.5, cmap=CMAP_IC,colorbar=False)
#plt.savefig(f'../results/Figures/IC_surf_r_lateral.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(IC_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map = col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_right'], bg_on_data = False,
darkness=0.5,cmap=CMAP_IC,colorbar=False)
#plt.savefig(f'../results/Figures/IC_surf_r_medial.png', dpi=300, bbox_inches='tight')
plt.show()
cmap = plt.get_cmap(CMAP_IC)
norm = plt.Normalize(vmin=V_min, vmax=V_max)
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
cb = colorbar.ColorbarBase(ax, cmap=cmap, norm=norm, orientation='horizontal')
cb.ax.tick_params(labelsize=18)
#plt.savefig('../results/Figures/IC_colorbar_horizontal.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#---------------------- regions with strongest IC values -----------------------------------
#--------------------------------------------------------------------------------------------
#plot_node_surf(DATA_avg, MIconn_allsub, 'ICallsub_w_avg', CMAP_IC,0.1, limb_ind)
# %% [markdown]
# ## Fig 2b: DC vs CMRglc
# %%
deg_color = '#9CB9E8'
COLOR = deg_color
corr_deg_w_E_allsub = np.zeros((sub_size))
p_deg_w_E_allsub = np.zeros((sub_size))
for j in range(sub_size):
corr_deg_w_E_allsub[j], p_deg_w_E_allsub[j] = pcor(degallsub_w_pr[:,j], edata_medianallsub[:, j])
corr_deg_w_E, p_deg_w_E = pcor(DATA_avg_pr.degallsub_w_avg , DATA_avg_pr.pet_avg)
#--------------------------------------------------------------------------------------------
#---------------------- spatial autocorrelation (SAC) -----------------------------------
#--------------------------------------------------------------------------------------------
### BRAI SMASH:
# niter = 1000
# test_stat_DC ,surrogate_brainmap_corrs_DC, sa_corrected_p_value_DC, spatially_naive_p_value_DC = Spatial_AC(DATA_avg_pr , "deg_w", LIMB,niter)
# sac1 = '#0199DD'
# sac2 = '#8CCF83'
# sac3 = '#FEAD01'
# fig, ax = plt.subplots(figsize=(3, 7))
# plt.grid(True)
# g = sns.kdeplot(surrogate_brainmap_corrs_DC, color= "#E0E0E0", ax=ax, fill= True, linewidth=2.5 )
# ax.axvline(corr_deg_w_E, 0, 0.9, color= COLOR, linestyle='dashed', lw=3)
# ax.set_xticks(np.arange(-1, 1.1, 0.5))
# ax.set_ylim(0, 2)
# ax.set_xlim(-0.8, 0.8)
# #ax.text(0.5, -0.1, "Pearson correlation\nwith IE map", ha='center', va='top', transform=ax.transAxes)
# spine_width = 1.6
# spine_color = 'gray'
# for spine in ax.spines.values():
# spine.set_linewidth(spine_width)
# spine.set_color(spine_color)
# plt.tick_params(axis='both', which='major', labelsize= FONT)
# plt.xlabel('Correlation', fontsize = FONT)
# plt.ylabel('Density', fontsize = FONT)
# plt.savefig(f'../results/Figures/CMRglc_DC_hist.png', dpi=300, bbox_inches='tight')
# plt.show()
# lower_bound = np.percentile(surrogate_brainmap_corrs_DC, 5)
# upper_bound = np.percentile(surrogate_brainmap_corrs_DC, 95)
# print(f"95% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
#--------------------------------------------------------------------------------------------
#---------------------- spin test -----------------------------------
#--------------------------------------------------------------------------------------------
N_ROT = 10000
RNG_SEED = 42
session = 'AUF'
LIMB = 'with'
data_avg_path_pr = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
DATA_avg_spin = pd.read_csv(data_avg_path_pr)
cmrglc = DATA_avg_spin['pet_avg'].to_numpy(dtype=float).ravel()
deg = DATA_avg_spin['degallsub_w_avg'].to_numpy(dtype=float).ravel()
rho_deg, p_spin_deg, null_deg = run_spin_test(deg, cmrglc, N_ROT,RNG_SEED)
FONT = 22
null_r = np.asarray(null_deg, dtype=float)
null_r = null_r[np.isfinite(null_r)]
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
sns.kdeplot(
null_r,
color="#E0E0E0",
ax=ax,
fill=True,
linewidth=2.5
)
ax.axvline(rho_deg, 0, 0.9, color=COLOR, linestyle='dashed', lw=3)
ax.set_xticks(np.arange(-1, 1.1, 0.5))
ax.set_ylim(0, 4)
ax.set_xlim(-0.8, 0.8)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.tick_params(axis='both', which='major', labelsize=FONT)
plt.xlabel('Correlation', fontsize=FONT)
plt.ylabel('Density', fontsize=FONT)
#plt.savefig('../results/Figures/CMRglc_DC_hist.png', dpi=300, bbox_inches='tight')
plt.show()
lower_bound = np.percentile(null_r, 5)
upper_bound = np.percentile(null_r, 95)
print(f"90% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
#--------------------------------------------------------------------------------------------
#------------------ scatter deg_avg , CMRglc_avg / group analysis ----------------------------
#--------------------------------------------------------------------------------------------
plt.figure(figsize=(6,4))
g = sns.jointplot(data=DATA_avg_pr, x='pet_avg', y='degallsub_w_avg', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
p_deg_w_E_smash = p_spin_deg
if p_deg_w_E_smash < 0.0001:
title = "r = {:.2f} , p-spin< {:.3f}".format(float(corr_deg_w_E), float(0.001))
else:
title = "r = {:.2f} , p-spin = {:.2f}".format(float(corr_deg_w_E), float(p_deg_w_E_smash))
#title = "corr = {:.2f} , pvalue = {:.0e}".format(float(corr_deg_w_E), float(p_deg_w_E))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'Degree', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Target Energy (CMRglc)', fontsize = FONT)
plt.ylabel('DC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig(f'../results/Figures/CMRglc_DC_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
# #--------------------------------------------------------------------------------------------
# #------------------ Box plot individual correlation between IC and E -----------------------
# #--------------------------------------------------------------------------------------------
data = pd.DataFrame({'corr': corr_deg_w_E_allsub, 'pval': p_deg_w_E_allsub})
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
violin = sns.violinplot(data=data, y='corr', ax=ax, color = COLOR)
violin.collections[0].set_facecolor('none')
sns.stripplot(data=data[data['pval'] < 0.05], y='corr', color=COLOR , size = 8, ax=ax)
sns.stripplot(data=data[data['pval'] >= 0.05], y='corr', color='red' , size = 8, ax=ax)
ax.set_ylabel('Correlation DC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.1,0.72])
plt.tick_params(axis='both', which='major', labelsize= FONT)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
x_position = 0.25
y_position = 0.92
ax.scatter(x_position, y_position, s=80, color='red', transform=ax.transAxes)
ax.text(x_position + 0.08, y_position, 'p > 0.05',
color='gray', fontsize=FONT, va='center', ha='left', transform=ax.transAxes)
plt.savefig(f'../results/Figures/CMRglc_DC_box.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- Degree map on the brain surface -----------------------------------
#--------------------------------------------------------------------------------------------
deg = np.array(DATA_avg_pr.degallsub_w_avg).reshape(-1, 1)
deg_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[limb_ind] = False
ind = ind[mask]
deg_allrois[ind] = np.abs(deg)
fig = plt.figure(figsize=(5, 5),dpi=300)
fsaverage = datasets.fetch_surf_fsaverage()
V_min = deg.min()
V_max = deg.max()
#V_max = 34
col_fsa = parcel_to_surface(deg_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP_DC,colorbar=False)
plt.savefig(f'../results/Figures/DC_surf_l_lateral.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(deg_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data = False,
darkness=0.5,cmap=CMAP_DC,colorbar=False)
plt.savefig(f'../results/Figures/DC_surf_l_medial.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(deg_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map = col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_right'], bg_on_data=False,
darkness=0.5, cmap=CMAP_DC,colorbar=False)
plt.savefig(f'../results/Figures/DC_surf_r_lateral.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(deg_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map = col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_right'], bg_on_data = False,
darkness=0.5,cmap=CMAP_DC,colorbar=False)
plt.savefig(f'../results/Figures/DC_surf_r_medial.png', dpi=300, bbox_inches='tight')
plt.show()
cmap = plt.get_cmap(CMAP_DC)
norm = plt.Normalize(vmin=V_min, vmax=V_max)
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
cb = colorbar.ColorbarBase(ax, cmap=cmap, norm=norm, orientation='horizontal')
cb.ax.tick_params(labelsize=18)
plt.savefig('../results/Figures/DC_colorbar_horizontal.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#---------------------- regions with strongest IC values -----------------------------------
#--------------------------------------------------------------------------------------------
#plot_node_surf(DATA_avg_pr, Pearconn_allsub, 'degallsub_w_avg', CMAP_DC,0.3, limb_ind)
# %% [markdown]
# ## Fig 3a: MwC and DC across networks
# %%
#--------------------------------------------------------------------------------------------
#--------------------------- test across networks ----------------------------------------------
#--------------------------------------------------------------------------------------------
import scipy.stats as stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from matplotlib.cbook import boxplot_stats
# Group data by network
groups_ic = [DATA_avg[DATA_avg['yeo_7_nw'] == nw]['ICallsub_w_avg'].dropna() for nw in DATA_avg['yeo_7_nw'].unique()]
groups_dc = [DATA_avg_pr[DATA_avg_pr['yeo_7_nw'] == nw]['degallsub_w_avg'].dropna() for nw in DATA_avg_pr['yeo_7_nw'].unique()]
# Run ANOVA test
anova_mwc_result = stats.f_oneway(*groups_ic)
anova_dc_result = stats.f_oneway(*groups_dc)
print(f"ANOVA for MwC across networks: F = {anova_mwc_result.statistic:.2f}, p = {anova_mwc_result.pvalue:.4f}")
print(f"ANOVA for DC across networks: F = {anova_dc_result.statistic:.2f}, p = {anova_dc_result.pvalue:.4f}")
# Run Tukey's HSD test for MwC
tukey_mwc = pairwise_tukeyhsd(DATA_avg['ICallsub_w_avg'], DATA_avg['yeo_7_nw'])
tukey_df_mwc = pd.DataFrame(data=tukey_mwc._results_table.data[1:], columns=tukey_mwc._results_table.data[0])
sig_results_mwc = tukey_df_mwc[tukey_df_mwc['p-adj'] < 0.05]
# Run Tukey's HSD test for DC
tukey_dc = pairwise_tukeyhsd(DATA_avg_pr['degallsub_w_avg'], DATA_avg_pr['yeo_7_nw'])
tukey_df_dc = pd.DataFrame(data=tukey_dc._results_table.data[1:], columns=tukey_dc._results_table.data[0])
sig_results_dc = tukey_df_dc[tukey_df_dc['p-adj'] < 0.05]
#--------------------------------------------------------------------------------------------
#--------------------------- DC across networks ----------------------------------------------
#--------------------------------------------------------------------------------------------
# Calculate median values for each network
median_values = DATA_avg.groupby('yeo_7_nw')['ICallsub_w_avg'].median()
# Sort networks by the median values
ordered_networks = median_values.sort_values().index
# Generate ordered palettes based on the original colors
ordered_palette = [lighter_color_map[network] for network in ordered_networks]
strip_palette = [network_color_map[network] for network in ordered_networks]
# Plot with the correct colors and sorted order
plt.figure(figsize=(6, 4), dpi=300)
ax = plt.gca()
ax.set_axisbelow(True)
# Violin plot with lighter colors for fills
violin_parts = sns.violinplot(
data=DATA_avg, x='yeo_7_nw', y='ICallsub_w_avg',
inner=None, linewidth=1.5, palette=ordered_palette,
width=0.9, density_norm='width', ax=ax,
order=ordered_networks # Sort the networks based on median values
)
# Adjust violin outlines to use the correct lighter colors
for pc, network in zip(violin_parts.collections, ordered_networks):
mpath = np.array(pc.get_paths()[0].vertices)
mpath[:, 0] = np.clip(mpath[:, 0], np.median(mpath[:, 0]), np.max(mpath[:, 0]))
pc.set_paths([mpath])
pc.set_edgecolor(lighter_color_map[network]) # Use the correct lighter color for outlines
# Strip plot with the correct network colors
strip_plot = sns.stripplot(
data=DATA_avg, x='yeo_7_nw', y='ICallsub_w_avg',
jitter=0.1, marker='o', alpha=0.7,
palette=strip_palette, ax=ax, size=3, order=ordered_networks
)
# Offset points in the strip plot (if required)
for patch in strip_plot.collections:
x_shift = -0.15
patch.set_offsets(np.c_[patch.get_offsets()[:, 0] + x_shift, patch.get_offsets()[:, 1]])
# Box plot with updated order but original style
sns.boxplot(
data=DATA_avg, x='yeo_7_nw', y='ICallsub_w_avg',
width=0.07, fliersize=0, ax=ax,
boxprops={'facecolor': 'none', 'edgecolor': '#595B61', 'linewidth': 1},
whiskerprops={'color': '#595B61', 'linewidth': 1},
capprops={'color': '#595B61', 'linewidth': 1},
medianprops={'color': '#595B61', 'linewidth': 1},
order=ordered_networks
)
xvals = list(ordered_networks)
networks = sorted(DATA_avg['yeo_7_nw'].unique())
networks = sorted(DATA_avg['yeo_7_nw'].unique())
box_top = {}
for net in networks:
vals = DATA_avg.loc[DATA_avg['yeo_7_nw'] == net, 'ICallsub_w_avg']
stats = boxplot_stats(vals, whis=1.5)[0]
whisker = stats["whishi"]
actual_max = vals.max()
# Use whichever is higher
box_top[net] = max(whisker, actual_max)
if net == 'Default':
box_top[net] =box_top[net] + 0.8
star_counts = {net: 0 for net in ordered_networks}
star_offset = 0.8
star_spacing = 0.3
for idx, row in sig_results_mwc.iterrows():
group1 = row['group1']
group2 = row['group2']
p_value = row['p-adj']
star_text = get_star(p_value)
x_coord1 = xvals.index(group1)
y_coord1 = box_top[group1] + star_offset + star_counts[group1] * star_spacing
ax.text(x_coord1, y_coord1, star_text,
ha='center', va='bottom', fontsize=8, color=network_color_map[group2])
star_counts[group1] += 1
x_coord2 = xvals.index(group2)
y_coord2 = box_top[group2] + star_offset + star_counts[group2] * star_spacing
ax.text(x_coord2, y_coord2, star_text,
ha='center', va='bottom', fontsize=8, color=network_color_map[group1])
star_counts[group2] += 1
for spine in ax.spines.values():
spine.set_edgecolor('gray')
plt.grid(True, axis='both', which='major', linestyle='--', linewidth=0.5, color='gray')
ax.set_xticklabels([])
plt.xlabel('')
plt.ylabel('MwC')
plt.ylim([24,39])
plt.tight_layout()
plt.savefig(f'../results/Figures/network_IC_psmash.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- DC across networks ----------------------------------------------
#--------------------------------------------------------------------------------------------
# Calculate median values for each network
median_values = DATA_avg_pr.groupby('yeo_7_nw')['degallsub_w_avg'].median()
# Sort networks by the median values
ordered_networks = median_values.sort_values().index
# Generate ordered palettes based on the original colors
ordered_palette = [lighter_color_map[network] for network in ordered_networks]
strip_palette = [network_color_map[network] for network in ordered_networks]
# Plot with the correct colors and sorted order
plt.figure(figsize=(6, 4), dpi=300)
ax = plt.gca()
# Violin plot with lighter colors for fills
violin_parts = sns.violinplot(
data=DATA_avg_pr, x='yeo_7_nw', y='degallsub_w_avg',
inner=None, linewidth=1.5, palette=ordered_palette,
width=0.9, density_norm='width', ax=ax,
order=ordered_networks # Sort the networks based on median values
)
# Adjust violin outlines to use the correct lighter colors
for pc, network in zip(violin_parts.collections, ordered_networks):
mpath = np.array(pc.get_paths()[0].vertices)
mpath[:, 0] = np.clip(mpath[:, 0], np.median(mpath[:, 0]), np.max(mpath[:, 0]))
pc.set_paths([mpath])
pc.set_edgecolor(lighter_color_map[network]) # Use the correct lighter color for outlines
# Strip plot with the correct network colors
strip_plot = sns.stripplot(
data=DATA_avg_pr, x='yeo_7_nw', y='degallsub_w_avg',
jitter=0.1, marker='o', alpha=0.7,
palette=strip_palette, ax=ax, size=3, order=ordered_networks
)
# Offset points in the strip plot (if required)
for patch in strip_plot.collections:
x_shift = -0.15
patch.set_offsets(np.c_[patch.get_offsets()[:, 0] + x_shift, patch.get_offsets()[:, 1]])
# Box plot with updated order but original style
sns.boxplot(
data=DATA_avg_pr, x='yeo_7_nw', y='degallsub_w_avg',
width=0.07, fliersize=0, ax=ax,
boxprops={'facecolor': 'none', 'edgecolor': '#595B61', 'linewidth': 1},
whiskerprops={'color': '#595B61', 'linewidth': 1},
capprops={'color': '#595B61', 'linewidth': 1},
medianprops={'color': '#595B61', 'linewidth': 1},
order=ordered_networks
)
xvals = list(ordered_networks)
networks = sorted(DATA_avg_pr['yeo_7_nw'].unique())
networks = sorted(DATA_avg_pr['yeo_7_nw'].unique())
box_top = {}
for net in networks:
print(net)
vals = DATA_avg_pr.loc[DATA_avg_pr['yeo_7_nw'] == net, 'degallsub_w_avg']
stats = boxplot_stats(vals, whis=1.5)[0]
whisker = stats["whishi"]
actual_max = vals.max()
# Use whichever is higher
box_top[net] = max(whisker, actual_max) +2
if net == 'Vis':
box_top[net] =box_top[net] + 1
star_counts = {net: 0 for net in ordered_networks}
star_offset = 0.8
star_spacing = 0.3
for idx, row in sig_results_dc.iterrows():
group1 = row['group1']
group2 = row['group2']
p_value = row['p-adj']
star_text = get_star(p_value)
x_coord1 = xvals.index(group1)
y_coord1 = box_top[group1] + star_offset + star_counts[group1] * star_spacing
ax.text(x_coord1, y_coord1, star_text,
ha='center', va='bottom', fontsize= 10, color=network_color_map[group2])
star_counts[group1] += 1
x_coord2 = xvals.index(group2)
y_coord2 = box_top[group2] + star_offset + star_counts[group2] * star_spacing
ax.text(x_coord2, y_coord2, star_text,
ha='center', va='bottom', fontsize=10, color=network_color_map[group1])
star_counts[group2] += 1
for spine in ax.spines.values():
spine.set_edgecolor('gray')
# Final plot adjustments
ax.set_xticklabels([])
plt.grid(True, linestyle='--', linewidth=0.5, color='gray')
plt.xlabel('')
plt.ylabel('DC')
plt.tight_layout()
plt.savefig(f'../results/Figures/network_DC_psmash.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig 3b: Diversity
# %%
#--------------------------------------------------------------------------------------------
#--------------------------- MwC diversity ----------------------------------------------
#--------------------------------------------------------------------------------------------
from scipy.stats import spearmanr
thr_sc = 0.1
network_colors = ['#9A199A', '#45B7F0','#58B73B','#DE4B52','#E237FF','#F29A36','#FEFFD3']
Conn_Mat_IC = MIconn_allsub.copy()
Conn_Mat_IC[Conn_Mat_IC < 0.001] = 0
FCmat = Conn_Mat_IC
SCmat = SCconn_allsub.copy()
SCmat_th = np.zeros((SCmat.shape[0], SCmat.shape[1], sub_size))
for i in range(sub_size):
SCmat_th[:, :, i] = threshold_proportional(SCmat[:, :, i], thr_sc)
SC_mask = SCmat_th.copy()
SC_mask[SC_mask > 0] = 1
FCmat_SC = np.multiply(FCmat , SC_mask)
net_assign_allsub_IC = np.zeros((nrois_rem,num_net, sub_size))
for i in range(sub_size):
FC = FCmat_SC[:,:,i]
FC[FC > 0] = 1
#FC[FC < 0] = 1
FC_net = np.multiply(FC, new_net_num)
for r in range(nrois_rem):
for s in range(num_net):
net_assign_allsub_IC[r,s,i] = np.nansum(FC_net[r,:] == s+1)
div_allsub_IC = np.zeros((nrois_rem, sub_size))
for i in range(sub_size):
net_assign = net_assign_allsub_IC[:,:,i]
conn_sum = np.nansum(net_assign, axis=1, keepdims=True)
probabilities = net_assign / conn_sum
div_allsub_IC[:,i] = 1 - np.nansum(probabilities ** 2 , axis = 1)
div_avg_IC = np.nanmean(div_allsub_IC , axis = 1)
#
# --------- finding p-smash --------------------------------------
# ic_array = DATA_avg['ICallsub_w_avg'].to_numpy()
# niter = 1000
# data = pd.DataFrame({"x":div_avg_IC, "y":ic_array})
# test_stat_ic ,surrogate_brainmap_corrs_ic, sa_corrected_p_value_ic, spatially_naive_p_value_ic = Spatial_AC(data , "div", LIMB,niter)
# sac1 = '#0199DD'
# sac2 = '#8CCF83'
# sac3 = '#FEAD01'
# fig, ax = plt.subplots(figsize=(3, 7))
# plt.grid(True)
# g = sns.kdeplot(surrogate_brainmap_corrs_ic, color= "#E0E0E0", ax=ax, fill= True, linewidth=2.5 )
# ax.axvline(test_stat_ic, 0, 0.96, color= COLOR, linestyle='dashed', lw=3)
# ax.set_xticks(np.arange(-1, 1.1, 0.5))
# ax.set_ylim(0, 1.8)
# ax.set_xlim(-0.8, 0.8)
# spine_width = 1.6
# spine_color = 'gray'
# for spine in ax.spines.values():
# spine.set_linewidth(spine_width)
# spine.set_color(spine_color)
# plt.tick_params(axis='both', which='major', labelsize= FONT)
# plt.xlabel('Correlation', fontsize = FONT)
# plt.ylabel('Density', fontsize = FONT)
# plt.show()
# lower_bound = np.percentile(surrogate_brainmap_corrs_ic, 5)
# upper_bound = np.percentile(surrogate_brainmap_corrs_ic, 95)
#print(f"95% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
# ---------------------------------------------------
plt.figure(figsize=(5, 5), dpi = 300)
cor , p = pcor(div_avg_IC , DATA_avg['ICallsub_w_avg'])
p_IC_w_E_smash = p
data = pd.DataFrame({'ent_avg': div_avg_IC, 'ICallsub_w_avg': DATA_avg['ICallsub_w_avg'], 'Network': net_label["yeo_7_nw"]})
grid = sns.FacetGrid(data, hue="Network", palette=network_colors, height=6, aspect=5/4)
grid.map(plt.scatter, 'ent_avg', 'ICallsub_w_avg', alpha=0.9 , s=85)
sns.regplot(data=data, x='ent_avg', y='ICallsub_w_avg', scatter=False, ax=grid.ax, color='gray')
plt.grid(True)
if p_IC_w_E_smash < 0.001:
title = "r = {:.2f} , p < {:.3f}".format(float(cor), float(0.001))
else:
title = "r = {:.2f} , p = {:.2f}".format(float(cor), float(p_IC_w_E_smash))
#plt.title(title , fontsize=FONT)
grid.ax.set_title(title, fontsize=FONT)
plt.xlabel('Simpson Diversity Index', fontsize=FONT)
plt.ylabel('MwC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize=FONT)
for _, spine in grid.ax.spines.items():
spine.set_visible(True)
spine.set_linewidth(1)
spine.set_edgecolor('gray')
plt.savefig(f'../results/Figures/diversity_IC_psmash.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- DC diversity ----------------------------------------------
#--------------------------------------------------------------------------------------------
thr_FC = 0.2
thr_SC = 0.1
Conn_Mat_pr = Pearconn_allsub.copy()
Conn_Mat_pr[Conn_Mat_pr < 0] = 0
Conn_Mat_th = np.zeros((Conn_Mat_pr.shape[0], Conn_Mat_pr.shape[1], sub_size))
for i in range(sub_size):
Conn_Mat_th[:,:,i] = threshold_proportional(Conn_Mat_pr[:, :, i], thr_FC)
SCmat = SCconn_allsub.copy()
SCmat_th = np.zeros((SCmat.shape[0], SCmat.shape[1], sub_size))
for i in range(sub_size):
SCmat_th[:, :, i] = threshold_proportional(SCmat[:, :, i], thr_SC)
SC_mask = SCmat_th.copy()
SC_mask[SC_mask > 0] = 1
FCmat_SC_pr = np.multiply(Conn_Mat_th , SC_mask)
degallsub = np.zeros((nrois_rem, sub_size))
for i in range(sub_size):
dd = FCmat_SC_pr [:,:,i]
degallsub[:,i] = np.nansum(dd, axis=1)
net_assign_allsub_pr = np.zeros((nrois_rem,num_net, sub_size))
for i in range(sub_size):
FC = np.abs(FCmat_SC_pr[:,:,i])
FC[FC > 0] = 1
FC_net = np.multiply(FC, new_net_num)
for r in range(nrois_rem):
for s in range(num_net):
net_assign_allsub_pr[r,s,i] = np.nansum(FC_net[r,:] == s+1)
div_allsub_DC = np.zeros((nrois_rem, sub_size))
for i in range(sub_size):
net_assign = net_assign_allsub_pr[:,:,i]
conn_sum = np.nansum(net_assign, axis=1, keepdims=True)
probabilities = net_assign / conn_sum
div_allsub_DC[:,i] = 1 - np.nansum(probabilities ** 2 , axis = 1)
div_avg_DC = np.nanmean(div_allsub_DC , axis = 1)
# --------- finding p-smash -----------------------------------------
deg_array = DATA_avg_pr['degallsub_w_avg'].to_numpy()
cor , p = pcor(div_avg_DC , deg_array)
# niter = 10
# data = pd.DataFrame({"x":div_avg_DC, "y":deg_array})
# sac1 = '#0199DD'
# sac2 = '#8CCF83'
# sac3 = '#FEAD01'
# fig, ax = plt.subplots(figsize=(3, 7))
# plt.grid(True)
# g = sns.kdeplot(surrogate_brainmap_corrs_dc, color= "#E0E0E0", ax=ax, fill= True, linewidth=2.5 )
# ax.axvline(test_stat_dc, 0, 0.96, color= COLOR, linestyle='dashed', lw=3)
# ax.set_xticks(np.arange(-1, 1.1, 0.5))
# ax.set_ylim(0, 1.8)
# ax.set_xlim(-0.8, 0.8)
# spine_width = 1.6
# spine_color = 'gray'
# for spine in ax.spines.values():
# spine.set_linewidth(spine_width)
# spine.set_color(spine_color)
# plt.tick_params(axis='both', which='major', labelsize= FONT)
# plt.xlabel('Correlation', fontsize = FONT)
# plt.ylabel('Density', fontsize = FONT)
# plt.show()
# lower_bound = np.percentile(surrogate_brainmap_corrs_dc, 5)
# upper_bound = np.percentile(surrogate_brainmap_corrs_dc, 95)
# print(f"95% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
# ----------------------------------------------------------------
plt.figure(figsize=(5, 5), dpi = 300)
data = pd.DataFrame({'ent_avg': div_avg_DC, 'degallsub_w_avg': deg_array, 'Network': net_label["yeo_7_nw"]})
grid = sns.FacetGrid(data, hue="Network", palette=network_colors, height=6, aspect=5/4)
grid.map(plt.scatter, 'ent_avg', 'degallsub_w_avg', alpha=0.9 , s=85)
sns.regplot(data=data, x='ent_avg', y='degallsub_w_avg', scatter=False, ax=grid.ax, color='gray')
p_DC_w_E_smash = p
if p_DC_w_E_smash < 0.001:
title = "r = {:.2f} , p < {:.3f}".format(float(cor), float(0.001))
else:
title = "r = {:.2f} , p = {:.2f}".format(float(cor), float(p_IC_w_E_smash))
plt.grid(True)
#plt.title(title, fontsize=FONT)grid.ax.set_title(title, fontsize=FONT)
grid.ax.set_title(title, fontsize=FONT)
plt.xlabel('Simpson Diversity Index', fontsize=FONT)
plt.ylabel('DC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize=FONT)
for _, spine in grid.ax.spines.items():
spine.set_visible(True)
spine.set_linewidth(1)
spine.set_edgecolor('gray')
plt.savefig(f'../results/Figures/diversity_DC_psmash.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig 3c: Neurosynth
# %%
alltasks_rois = pd.read_csv('../data/external/Neurosynth_roi.csv', header=None)
alltasks_rois_rem = np.delete(alltasks_rois, limb_ind, axis =0 )
labels = ['face', 'verbal semantics', 'cued attention', 'working memory','autobiographical memory', 'reading', 'inhibition', 'motor',
'visual perception', 'numerical cognition', 'reward', 'visual attention','multisensory', 'visuospatial','eye movements', 'action',
'auditory', 'pain', 'language', 'declarative memory','visual semantics', 'emotion', 'cognitive control', 'social cognition']
base_radius = 0.25
correlation_values = np.linspace(0, 0.4, 5)
max_correlation = max(correlation_values)
#--------------------------------------------------------------------------------------------
#------------------ Neurosynth and IC map -------------------------------
#-------------------------------------------------------------------------------------------
CMAP = 'viridis'
#data preperation
data = DATA_avg.ICallsub_w_avg
IC = np.array((data)).reshape(-1, 1)
IC_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[limb_ind] = False
ind = ind[mask]
IC_allrois[ind] = (np.abs(IC))
cor =np.zeros((len(labels)))
p =np.zeros((len(labels)))
for i, label in enumerate(labels):
cor[i],p[i] = pcor(alltasks_rois_rem[:,i],data )
data = pd.DataFrame({'corr':cor, 'Pval':p, 'task':labels})
labels = data['task']
stats = data['corr']
Pval = data['Pval']
df = pd.DataFrame({'label': labels, 'corr': stats, 'Pval':Pval})
df_sorted = df.sort_values('corr', ascending=False).reset_index(drop=True)
sorted_labels = df_sorted['label'].tolist()
sorted_stats = df_sorted['corr'].tolist()
sorted_p = df_sorted['Pval'].tolist()
V_max = IC.max()
V_max = 34
V_min = IC.min()
fig = plt.figure(figsize=(6, 6))
ax_brain = fig.add_axes([0.38, 0.38, 0.24, 0.24], projection='3d')
fsaverage = datasets.fetch_surf_fsaverage()
col_fsa = parcel_to_surface(IC_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP_IC, axes=ax_brain, figure=fig, colorbar=False)
ax_brain.set_axis_off()
ax = fig.add_subplot(111, polar=True, label="PolarPlot", frame_on=False)
# Generate angles for each sector
angles = np.linspace(0, 2 * np.pi, len(sorted_labels) + 1, endpoint=True).tolist()
# Compute sector centers by averaging adjacent angles
sector_centers = [(angles[i] + angles[i + 1]) / 2 for i in range(len(angles) - 1)]
for angle in angles:
ax.plot([angle, angle], [base_radius, base_radius + max_correlation + 0.05],
linestyle='dashed', color='gray', linewidth=0.5)
# Bar width
width = (2 * np.pi / len(sorted_labels))
# Base radius for bars
base_radius = 0.25
# Plot bars
bars = ax.bar(angles[:-1], np.abs(sorted_stats), width=width, bottom=base_radius,
align='edge', edgecolor='white')
# Set colors based on correlation values
for bar, stat, pval in zip(bars, sorted_stats, sorted_p):
bar.set_facecolor('#F5A614' if (stat > 0 and pval < 0.05) else ('#35B0F7' if (stat < 0 and pval < 0.05) else 'gray'))
bar.set_alpha(0.9)
# Correlation ring values
correlation_values = np.linspace(0, 0.4, 5)
max_correlation = max(correlation_values)
for r in correlation_values:
ax.plot(np.linspace(0, 2 * np.pi, 100), [base_radius + r] * 100, '--', color='gray', linewidth=0.5)
ax.text(0, base_radius + r , f'{r:.1f}', horizontalalignment='center', verticalalignment='bottom', fontsize = 10)
# Label placement adjustments
dash_line_length = base_radius + max_correlation + 0.05
label_distance = dash_line_length + 0.02
for angle, label in zip(sector_centers, sorted_labels):
if label == "autobiographical memory":
label = "autobiog memory"
if len(label) > 10:
label = label.replace(' ', '\n', 1)
# alignment based on angle
c = np.cos(angle) # + on right half, – on left half, ~0 at top/bottom
eps = 0.15 # how wide to treat as “vertical” (tune 0.10–0.20)
if abs(c) < eps: # near top/bottom → center
ha = 'center'
elif c > 0: # right half → left-align outward
ha = 'left'
else: # left half → right-align outward
ha = 'right'
ax.text(angle, label_distance, label,
ha=ha, va='center', fontsize=10, color='black', clip_on=False)
ax.set_xticks([])
ax.set_yticklabels([])
ax.grid(False)
ax.spines['polar'].set_visible(False)
plt.savefig(f'../results/Figures/neurosynth_IC.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------ Neurosynth and DC map -------------------------------
#--------------------------------------------------------------------------------------------
# cool_cmap = cm.get_cmap('cool')
# cool_colors = cool_cmap(np.linspace(0.15, 0.95, 256))
# CMAP = LinearSegmentedColormap.from_list("adjusted_cool", cool_colors)
#data preperation
data = DATA_avg_pr.degallsub_w_avg
DC = np.array((data)).reshape(-1, 1)
DC_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[limb_ind] = False
ind = ind[mask]
DC_allrois[ind] = (np.abs(DC))
cor =np.zeros((len(labels)))
p =np.zeros((len(labels)))
for i, label in enumerate(labels):
cor[i],p[i] = pcor(alltasks_rois_rem[:,i],data )
data = pd.DataFrame({'corr':cor, 'Pval':p, 'task':labels})
labels = data['task']
stats = data['corr']
Pval = data['Pval']
df = pd.DataFrame({'label': labels, 'corr': stats, 'Pval':Pval})
df_sorted = df.sort_values('corr', ascending=False).reset_index(drop=True)
sorted_labels = df_sorted['label'].tolist()
sorted_stats = df_sorted['corr'].tolist()
sorted_p = df_sorted['Pval'].tolist()
V_max = DC.max()
V_min = DC.min()
fig = plt.figure(figsize=(6, 6))
ax_brain = fig.add_axes([0.38, 0.38, 0.24, 0.24], projection='3d')
fsaverage = datasets.fetch_surf_fsaverage()
col_fsa = parcel_to_surface(DC_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP_DC, axes=ax_brain, figure=fig, colorbar=False)
ax_brain.set_axis_off()
ax = fig.add_subplot(111, polar=True, label="PolarPlot", frame_on=False)
# Generate angles for each sector
angles = np.linspace(0, 2 * np.pi, len(sorted_labels) + 1, endpoint=True).tolist()
# Compute sector centers by averaging adjacent angles
sector_centers = [(angles[i] + angles[i + 1]) / 2 for i in range(len(angles) - 1)]
for angle in angles:
ax.plot([angle, angle], [base_radius, base_radius + max_correlation + 0.05],
linestyle='dashed', color='gray', linewidth=0.5)
# Bar width
width = (2 * np.pi / len(sorted_labels))
# Base radius for bars
# Plot bars
bars = ax.bar(angles[:-1], np.abs(sorted_stats), width=width, bottom=base_radius,
align='edge', edgecolor='white')
# Set colors based on correlation values
for bar, stat, pval in zip(bars, sorted_stats, sorted_p):
bar.set_facecolor('#F5A614' if (stat > 0 and pval < 0.05) else ('#35B0F7' if (stat < 0 and pval < 0.05) else 'gray'))
bar.set_alpha(0.9)
# Correlation ring values
correlation_values = np.linspace(0, 0.4, 5)
max_correlation = max(correlation_values)
for r in correlation_values:
ax.plot(np.linspace(0, 2 * np.pi, 100), [base_radius + r] * 100, '--', color='gray', linewidth=0.5)
ax.text(0, base_radius + r, f'{r:.1f}', horizontalalignment='center', verticalalignment='bottom',
fontsize=10)
# Label placement adjustments
dash_line_length = base_radius + max_correlation + 0.05
label_distance = dash_line_length + 0.02
for angle, label in zip(sector_centers, sorted_labels):
if label == "multisensory":
label = "multi-\nsensory"
if label == "autobiographical memory":
label = "autobiog memory"
if len(label) > 9:
label = label.replace(' ', '\n', 1)
# alignment based on angle
c = np.cos(angle) # + on right half, – on left half, ~0 at top/bottom
eps = 0.15 # how wide to treat as “vertical” (tune 0.10–0.20)
if abs(c) < eps: # near top/bottom → center
ha = 'center'
elif c > 0: # right half → left-align outward
ha = 'left'
else: # left half → right-align outward
ha = 'right'
ax.text(angle, label_distance, label,
ha=ha, va='center', fontsize=10, color='black', clip_on=False)
ax.set_xticks([])
ax.set_yticklabels([])
ax.grid(False)
ax.spines['polar'].set_visible(False)
plt.savefig(f'../results/Figures/neurosynth_DC.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig 4a: PLS
# %%
genes = pd.read_csv("../data/external/glasser_expression.csv" )
gene_names = genes.columns[0:]
CMAP = 'plasma'
nregs_lh = 167
limb_ind_l = np.array([88, 90, 92, 93, 110, 118, 120, 122, 131, 135, 165, 166, 172])-1
gene_exp_l = genes.to_numpy()[:180 , :]
gene_exp_l = np.delete(gene_exp_l , limb_ind_l , axis = 0)
rows_with_nan = np.any(np.isnan(gene_exp_l), axis=1)
rows_with_nan_ind = [i for i, val in enumerate(rows_with_nan) if val]
allgenes_expression = pd.read_csv("../data/external/glasser_expression.csv" )
#rows_with_nan_ind = []
data_avg = DATA_avg
nregs_lh = 167
IC_l = np.array(data_avg.ICallsub_w_avg[:nregs_lh])
IC_l = np.delete(IC_l , rows_with_nan_ind , axis = 0)
IC_l = np.array(IC_l).reshape(-1, 1)
SIZE = 50
#--------------------------------------------------------------------------------------------
#------------------------ IC map -----------------------------
#--------------------------------------------------------------------------------------------
CMAP = 'plasma'
gene_exp_l = genes.to_numpy()[:180 , :]
rows_with_nan = np.any(np.isnan(gene_exp_l), axis=1)
rows_with_nan_ind = [i for i, val in enumerate(rows_with_nan) if val]
data_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[180:] = False
mask[limb_ind_l] = False
mask[rows_with_nan_ind] = False
ind = ind[mask]
IC_l = zscore(IC_l)
data_allrois[ind] = IC_l
vmax1 = 2
vmin1 = -2
# vmax1 = IC_l.max()
# vmin1 = IC_l.min()
col_fsa = parcel_to_surface(data_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= vmin1, vmax = vmax1, view='lateral',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP,colorbar=False)
plt.savefig('../results/Figures/IC_l_lateral_GO.png', dpi=300, bbox_inches='tight')
plt.show()
col_fsa = parcel_to_surface(data_allrois,'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= vmin1, vmax = vmax1, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data = False,
darkness=0.5,cmap=CMAP,colorbar=False)
plt.savefig('../results/Figures/IC_l_medial_GO.png', dpi=300, bbox_inches='tight')
plt.show()
cmap = plt.get_cmap(CMAP)
norm = plt.Normalize(vmin=vmin1, vmax=vmax1)
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
cb = colorbar.ColorbarBase(ax, cmap=cmap, norm=norm, orientation='horizontal')
cb.ax.tick_params(labelsize=16)
plt.savefig('../results/Figures/IC_colorbar_GO.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------------- AHBA genes PLS1 scores -----------------------------------------
#--------------------------------------------------------------------------------------------
pls1_score = pd.read_csv('../data/external/PLS1_ROIscores.csv', header=None)
pls1_score = np.array(pls1_score).reshape(-1, 1)
pls1_score_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[180:] = False
mask[limb_ind_l] = False
mask[rows_with_nan_ind] = False
ind = ind[mask]
pls1_score = zscore(pls1_score)
pls1_score_allrois[ind] = pls1_score
# vmax1 = pls1_score.max()
# vmin1 = pls1_score.min()
vmax1 = 2
vmin1 = -2
fsaverage = datasets.fetch_surf_fsaverage()
col_fsa = parcel_to_surface(pls1_score_allrois, 'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map=col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= vmin1, vmax= vmax1, view='lateral',
bg_map=fsaverage['sulc_right'], bg_on_data= False,
darkness=0.5, cmap=CMAP, colorbar= False)
plt.savefig('../results/Figures/pls1_l_lateral_GO.png', dpi=300, bbox_inches='tight')
plt.show()
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map=col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= vmin1, vmax= vmax1, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP, colorbar= False)
plt.savefig('../results/Figures/pls1_l_medial_GO.png', dpi=300, bbox_inches='tight')
plt.show()
cmap = plt.get_cmap(CMAP)
norm = plt.Normalize(vmin=vmin1, vmax=vmax1)
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
cb = colorbar.ColorbarBase(ax, cmap=cmap, norm=norm, orientation='horizontal')
cb.ax.tick_params(labelsize=16)
plt.savefig('../results/Figures/IC_colorbar_GO.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------------- scatter plot IC and pls1 score ---------------------------------------
#--------------------------------------------------------------------------------------------
IC = IC_l.flatten()
pls1 = pls1_score.flatten()
corr_IC_w_pls1, p_IC_w_pls1 = pcor(IC, pls1 )
fig, ax = plt.subplots(figsize=(4, 4))
sns.regplot(
x=IC,
y=pls1,
scatter_kws={
'facecolors': "#E6E6E6",
'edgecolor': '#787B76' ,
'linewidths': 1,
's': 85
},
line_kws={
'color': '#000000',
'lw': 2
},
ci=None # This removes the confidence interval
)
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_w_pls1) , float(p_IC_w_pls1))
ax.set_title(title , fontdict={'fontsize': FONT})
plt.xlabel('Group level MwC' , fontsize = FONT)
plt.ylabel('PLS1 score' , fontsize = FONT)
plt.grid(True)
plt.tick_params(axis='both', which='major', labelsize= FONT) # Increase the tick label font size
ax = plt.gca() # Get the current Axes instance
# Set the spine linewidth
spine_width = 1.6 # Change this value to your preferred linewidth
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.savefig('../results/Figures/IC_pls1_GO.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig 4b: Gene enrichment
# %%
target_strings = ['axon', 'synap', 'example', 'dendritic', 'mito', 'signal' , 'transporter' , 'carbo'] # Add your specific strings here
FONT = 12
data = pd.read_csv(f'../data/external/metascape_result_summary.csv', delimiter=';')
pattern = '|'.join(target_strings)
data = data[data['Description'].str.contains(pattern, case=False, na=False)]
data['LogP'] = pd.to_numeric(data['Log10(P)'], errors='coerce')
data['-Log10(P)'] = -data['Log10(P)']
data['Description'] = data['Description'].apply(lambda x: x[0].upper() + x[1:] if pd.notnull(x) and len(x) > 0 else x)
p_values = data['LogP']
num_genes = data['Count']
descriptions = data['Description']
colors = ['#E4191D', '#367EB8', '#4DAF4A','#FFFC32','#F781BF','#999999','#8DD3C7','#FB8072','#D9D9D9'] # Add your desired colors here
target_strings = ['axon', 'synap', 'example', 'dendritic', 'mito', 'signal' , 'transporter' , 'carbo'] # Add your specific strings here
DPI = 300
min_count = data['Count'].min()
max_count = data['Count'].max()
sizes = (min_count*2.5, max_count*2.5)
start = min_count * 2.5
end = max_count * 2.5
step_size = (end - start) / 3
representative_sizes = [start, start + step_size, start + 2 * step_size, end]
plt.figure(figsize=(5, 6), dpi=DPI)
ax = plt.gca()
ax.set_axisbelow(True)
ax.grid(True)
g = sns.scatterplot(
x='-Log10(P)',
y='Description',
size='Count',
hue='Description',
sizes=(start, end),
palette=colors,
alpha=1,
data=data
)
# Remove seaborn’s legend
g.get_legend().remove()
# Wrap long descriptions onto multiple lines
y_labels = [textwrap.fill(desc, 30) for desc in data['Description']]
plt.yticks(range(len(y_labels)), y_labels, fontsize=FONT, fontweight='bold')
plt.xlabel('-log10(p)', fontsize=FONT+5)
plt.ylabel('', fontsize=FONT)
ax.set_axisbelow(True)
plt.xlim([0, 40])
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.6)
# --- Inset legend above the plotting rectangle ----------
legend_ax = ax.inset_axes([0, 1.03, 1, 0.08], transform=ax.transAxes)
legend_ax.axis('off')
legend_ax.set_zorder(10)
for size in representative_sizes:
legend_ax.scatter([], [], s=size, color='grey', alpha=0.6,
label=f'{int(size/2.5)}', zorder=11)
# Single legend call, centered in the inset, with smaller text
legend = legend_ax.legend(
title="Genes",
loc='center',
bbox_to_anchor=(0.5, 0.5),
ncol=len(representative_sizes),
columnspacing=0.8,
handletextpad=0.5,
frameon=False,
prop={'size': FONT-2}
)
legend.set_title("Genes", prop={'size': FONT-2})
legend._legend_box.align = "center"
# --- Save and show ----------
plt.savefig('../results/Figures/Gene_enrichment.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig 4c: synaptic density
# %%
mmp_in_mni_6 = f'../data/external/mmp_in_mni_6.nii.gz'
atlas_mni6 = nib.load(mmp_in_mni_6)
atlas_data_mni6 = atlas_mni6.get_fdata()
Bmax_mean_mni_img = nib.load('PATH to synaptic density atlas')
Bmax_mean_mni_data = Bmax_mean_mni_img.get_fdata()
region_labels = np.unique(atlas_data_mni6)
region_labels = region_labels[region_labels > 0]
thr = 0.001
n_regions = len(region_labels)
syn_values = np.zeros(n_regions)
for i, label in enumerate(region_labels):
region_mask = atlas_data_mni6 == label
region_vals_syn = Bmax_mean_mni_data[region_mask]
#syn_values[i] = np.median(region_vals_syn)
syn_values[i] = np.nanmedian(region_vals_syn[region_vals_syn > thr])
syn_values = np.delete(syn_values, limb_ind , axis = 0)
CMAP_IC = 'plasma'
# IC = np.array(DATA_avg.ICallsub_w_avg).reshape(-1, 1)
# V_min = 27
# V_max = 34
IC = np.array(syn_values).reshape(-1, 1)
V_min = 450
V_max = 700
IC_allrois = np.zeros((360,1))
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
mask[limb_ind] = False
ind = ind[mask]
IC_allrois[ind] = np.abs(IC)
fig = plt.figure(figsize=(5, 5),dpi=300)
fsaverage = datasets.fetch_surf_fsaverage()
col_fsa = parcel_to_surface(IC_allrois,'glasser_360_fsa5')
col_fsa[col_fsa == 0] = np.nan
col_fsa = np.where(np.isnan(col_fsa), np.nan, np.clip(col_fsa, V_min, V_max))
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
darkness=0.5, cmap=CMAP_IC,colorbar=False)
plt.savefig('../results/Figures/synaptic_den_l_lateral.png', dpi=300, bbox_inches='tight')
plt.show()
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map = col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data = False,
darkness=0.5,cmap=CMAP_IC,colorbar=False)
plt.savefig('../results/Figures/synaptic_den_l_medial.png', dpi=300, bbox_inches='tight')
plt.show()
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map = col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_right'], bg_on_data=False,
darkness=0.5, cmap=CMAP_IC,colorbar=False)
plt.savefig('../results/Figures/synaptic_den_r_lateral.png', dpi=300, bbox_inches='tight')
plt.show()
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map = col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_right'], bg_on_data = False,
darkness=0.5,cmap=CMAP_IC,colorbar=False)
plt.savefig('../results/Figures/synaptic_den_r_medial.png', dpi=300, bbox_inches='tight')
plt.show()
ticks = np.linspace(V_min, V_max, 6)
norm = plt.Normalize(vmin= V_min, vmax = V_max)
sm = plt.cm.ScalarMappable(cmap=CMAP, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, ticks=ticks, orientation='horizontal', pad=0.0001, shrink=0.6)
cmap = plt.get_cmap(CMAP_IC)
norm = plt.Normalize(vmin=V_min, vmax=V_max)
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
cb = colorbar.ColorbarBase(ax, cmap=cmap, norm=norm, orientation='horizontal')
cb.ax.tick_params(labelsize=16)
plt.savefig('../results/Figures/synaptic_den_colorbar.png', dpi=300, bbox_inches='tight')
plt.show()
#==================================================================
#==================================================================
ind = []
ind = syn_values>400
data = pd.DataFrame({'IC':DATA_avg.ICallsub_w_avg[ind],'syn':syn_values[ind],
'pet':DATA_avg.pet_avg[ind], 'deg':DATA_avg.degallsub_w_avg[ind], 'net':net_label['yeo_7_nw'][ind]})
metric = 'IC'
r, p = pearsonr(data['syn'], data[metric])
fig, ax = plt.subplots(figsize=(4,4))
FONT = 18
sns.regplot(
data = data,
x= 'syn',
y= 'IC',
scatter_kws={
'facecolors': "#E6E6E6",
'edgecolor': '#787B76' ,
'linewidths': 1,
's': 85
},
line_kws={
'color': '#000000',
'lw': 2
},
ci=None # This removes the confidence interval
)
title = "r = {:.2f} , p = {:.0e}".format(float(r) , float(p))
ax.set_title(title , fontdict={'fontsize': FONT})
plt.xlabel('Synaptic Density' , fontsize = FONT)
plt.ylabel('MwC' , fontsize = FONT)
plt.grid(True)
plt.tick_params(axis='both', which='major', labelsize= FONT) # Increase the tick label font size
ax = plt.gca() # Get the current Axes instance
# Set the spine linewidth
spine_width = 1.6 # Change this value to your preferred linewidth
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
plt.savefig('../results/Figures/IC_synaptic_den_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig 4d: KEGG and beta amyloid
# %%
data = pd.read_csv(f'../data/external/Enrichment_Analysis_Visualizer_data_KEGG_enrichr.csv')
print(data)
brain_terms = ['Oxidative phosphorylation','Pathways of neurodegeneration','Huntington disease',
'Parkinson disease','Prion disease','Thermogenesis','Alzheimer disease','Amyotrophic lateral sclerosis']
neurodegen_terms = ['Pathways of neurodegeneration','Huntington disease',
'Parkinson disease','Prion disease','Alzheimer disease','Amyotrophic lateral sclerosis']
subset_data = data.head(10)
subset_data = subset_data.iloc[::-1].reset_index(drop=True)
subset_data['NegLog10P'] = -np.log10(subset_data['p-value'])
plt.figure(figsize=(10, 7), dpi = 300)
bars = plt.barh(subset_data['term'], subset_data['NegLog10P'], color='#B3E2CD')
plt.yticks([])
gray = '#7A7777'
dark_gray = '#68686A'
black = 'black'
for bar, term, value in zip(bars, subset_data['term'], subset_data['p-value']):
if term in neurodegen_terms:
bar.set_edgecolor("#2E3834")
bar.set_linewidth(2) # adjust thickness as you like
else:
bar.set_edgecolor("none") # keep others without border
text_color = black if term in brain_terms else gray
fontweight = 550 if term in neurodegen_terms else 'normal'
Font = 21 if term in neurodegen_terms else 17
plt.text(bar.get_width() - (0.01 * max(subset_data['NegLog10P'])), bar.get_y() + bar.get_height() / 2,
term, va='center', ha='right', color=text_color, fontsize=Font)
plt.xlabel('-log10(p-value)' , fontsize = FONT+4)
plt.xticks(fontsize=FONT+4)
plt.tight_layout()
spine_width = 1.6 # Change this value to your preferred linewidth
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
#plt.savefig('../results/Figures/KEGG.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#---------------------- beta-amyloid -----------------------------------
#--------------------------------------------------------------------------------------------
limb_ind = np.array([88, 90, 92, 93, 110, 118, 120, 122, 131, 135, 165, 166, 172, 268, 270, 272, 273, 290, 298, 300, 302, 311, 314, 315, 346, 352])
roi_label = pd.read_csv("../data/external/hcp_mmp10_yeo7_modes.csv", delimiter=";")
roi_label["ROI_Name"] = (
roi_label["ROI_Name"]
.str.replace(r"^R_", "r", regex=True)
.str.replace(r"^L_", "l", regex=True)
.str.replace(r"_ROI$", "", regex=True)
.str.replace("-", ".", regex=False)
)
centrality ='MwC' # it can be DC or MCC
amyloid_status = 'Ab.pos' # it can be Ab.neg, Ab.pos
DX = 'MCI' #it can be MCI, CN, Dementia
remove_limbic = True
tracer = 'amyloid'
group_name = f"{DX}_{amyloid_status}"
group_data = pd.read_csv("../data/external/group_amyloid_mmp.csv")
amyloid = group_data[["ROI", group_name]].rename(columns={"ROI": "label", group_name: "Amyloid"})
roi = roi_label["ROI_Name"].reset_index(drop=True)
limb_labels = roi_label.loc[roi_label["ROI_Number"].isin(limb_ind + 1), "ROI_Name"]
if remove_limbic:
limb_labels = roi_label.loc[roi_label["ROI_Number"].isin(limb_ind + 1), "ROI_Name"]
amyloid = amyloid[~amyloid["label"].isin(limb_labels)].reset_index(drop=True)
roi = roi.drop(limb_ind, errors="ignore").reset_index(drop=True)
amyloid = amyloid.set_index("label").reindex(roi).reset_index()
IC = DATA_avg.ICallsub_w_avg if centrality == "MwC" else DATA_avg_pr.degallsub_w_avg
ad_IC = pd.DataFrame({
"IC": np.asarray(IC),
"Amyloid": amyloid["Amyloid"].values,
"label": roi.values,
"Network": net_label["yeo_7_nw"].reset_index(drop=True).values
})
corr_ic_ad, p_ic_ad = pcor(ad_IC['IC'], ad_IC['Amyloid'])
print("Correlation coefficient:", corr_ic_ad)
print("P-value:", p_ic_ad)
COLOR = "#8DD084"
plt.figure(figsize=(4,4))
plt.title = plt.gca().set_title
ax = sns.regplot(
x=ad_IC['Amyloid'],
y=ad_IC['IC'],
scatter_kws={
'facecolors': "#FFFFFF",
'edgecolor': '#43B180' ,
'linewidths': 2,
's': 85
},
line_kws={
'color': '#000000',
'lw': 1
},
ci=None # This removes the confidence interval
)
plt.title("r = {:.2f} , p-val = {:.0e}".format(float(corr_ic_ad), float(p_ic_ad)), fontsize=FONT)
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.0001", fontsize=FONT, va='baseline', y=0.8)
#g.fig.suptitle(f"corr = {corr_ic_ad:.2f}, Pval ={p_ic_ad:.5f}", fontsize=FONT, va='baseline', y=0.8)
#g.set_axis_labels('Amyloid', 'MCC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel(tracer + ' (' + amyloid_status + ', ' + DX + ')', fontsize=FONT)
plt.rcParams['axes.edgecolor'] = plt.rcParams['grid.color']
plt.ylabel(centrality, fontsize = FONT)
plt.savefig('../results/Figures/IC_bamyloid_MCI.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S1: hub details
# %%
plt.rcdefaults()
color_shades1 = ['#00915C','#519470','#748078']
color_shades2 = ['#E6AC00', '#b29049', '#7b7469']
color_shades3 = ['#0062F5', '#568FE3', '#8395AD']
data = pd.read_csv("../data/external/HCP-MMP1_UniqueRegionList.csv")
coord = data[["x-cog",'y-cog', 'z-cog']]
coord = coord.drop(limb_ind)
coord = coord.reset_index(drop = True)
roi_info = data[['cortex']]
roi_info= roi_info.drop(limb_ind)
roi_info = roi_info.reset_index(drop = True)
roi_label = data[['regionName']]
roi_label = roi_label.drop(limb_ind)
roi_label = roi_label.reset_index(drop = True)
roi_label1 = data[['Lobe']]
roi_label1 = roi_label1.drop(limb_ind)
roi_label1 = roi_label1.reset_index(drop = True)
roi_info["Network"] = net_label["yeo_7_nw"]
roi_info["Region"] = roi_label['regionName']
roi_info["Lobe"] = roi_label1['Lobe']
ICmap = DATA_avg.ICallsub_w_avg
wdegmap = DATA_avg_pr.degallsub_w_avg
ICmap = np.array(ICmap).reshape(-1, 1)
wdegmap = np.array(wdegmap).reshape(-1, 1)
ind = np.arange(360)
mask = np.ones(ind.shape, dtype=bool)
mask[limb_ind] = False
ind = ind[mask]
#--------------------------------------------------------------------------------------------
#---------------------- regions with strongest MwC values -----------------------------------
#--------------------------------------------------------------------------------------------
n_shades = 3
ic_cmap = LinearSegmentedColormap.from_list("ic_shades", color_shades1, N=n_shades)
#percentiles = [95, 90, 85, 80, 75]
percentiles = [95, 85, 75]
data_values = np.linspace(1, 5, len(percentiles))
data = np.zeros_like(ICmap)
for i, percentile in enumerate(percentiles):
threshold = np.percentile(ICmap, percentile)
mask = (ICmap >= threshold) & (data == 0)
data[mask] = i + 1
if percentile == 95:
data_ic = roi_info.loc[mask,:]
ind = np.arange(360)
mask = np.ones(ind.shape, dtype=bool)
mask[limb_ind] = False
ind = ind[mask]
data_allrois = np.zeros((360,1))
data_allrois[ind] = data
norm = Normalize(vmin=data_values.min(), vmax=data_values.max())
fsaverage = datasets.fetch_surf_fsaverage()
col_fsa = parcel_to_surface(data_allrois, 'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map=col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', view='medial', bg_map=fsaverage['sulc_left'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False,
title='MwC map')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map=col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', view='lateral', bg_map=fsaverage['sulc_left'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False)
legend_elements = [
Patch(facecolor=color_shades1[0], label='Top 5%'),
# Patch(facecolor=color_shades1[1], label='Top 10%'),
Patch(facecolor=color_shades1[1], label='Top 15%'),
# Patch(facecolor=color_shades1[3], label='Top 20%'),
Patch(facecolor=color_shades1[2], label='Top 25%')
]
plt.legend(handles=legend_elements, loc='upper center', bbox_to_anchor=(0.5, -0.05),
title='Percentile Ranges', ncol=len(legend_elements), frameon=False, fancybox=False, framealpha=0)
plt.show()
col_fsa = parcel_to_surface(data_allrois, 'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map=col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', view='medial', bg_map=fsaverage['sulc_right'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False,
title='MwC map')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map=col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', view='lateral', bg_map=fsaverage['sulc_right'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False)
legend_elements = [
Patch(facecolor=color_shades1[0], label='Top 5%'),
# Patch(facecolor=color_shades1[1], label='Top 10%'),
Patch(facecolor=color_shades1[1], label='Top 15%'),
# Patch(facecolor=color_shades1[3], label='Top 20%'),
Patch(facecolor=color_shades1[2], label='Top 25%')
]
plt.legend(handles=legend_elements, loc='upper center', bbox_to_anchor=(0.5, -0.05),
title='Percentile Ranges', ncol=len(legend_elements), frameon=False, fancybox=False, framealpha=0)
plt.show()
#--------------------------------------------------------------------------------------------
#---------------------- regions with strongest DC values -------------------------------
#--------------------------------------------------------------------------------------------
n_shades = 3
ic_cmap = LinearSegmentedColormap.from_list("ic_shades", color_shades3, N=n_shades)
percentiles = [95, 85, 75]
data_values = np.linspace(1, 5, len(percentiles))
data = np.zeros_like(wdegmap)
for i, percentile in enumerate(percentiles):
threshold = np.percentile(wdegmap, percentile)
mask = (wdegmap >= threshold) & (data == 0)
data[mask] = i + 1
if percentile == 95:
data_degree = roi_info.loc[mask,:]
#print( data_degree)
data_allrois = np.zeros((360,1))
data_allrois[ind] = data
norm = Normalize(vmin=data_values.min(), vmax=data_values.max())
fsaverage = datasets.fetch_surf_fsaverage()
col_fsa = parcel_to_surface(data_allrois, 'glasser_360_fsa5')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map=col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', view='medial', bg_map=fsaverage['sulc_left'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False,
title='DC map')
plotting.plot_surf_roi(fsaverage['infl_left'], roi_map=col_fsa[:int(col_fsa.shape[0]/2)],
hemi='left', view='lateral', bg_map=fsaverage['sulc_left'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False)
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map=col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', view='medial', bg_map=fsaverage['sulc_right'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False,
title='DC map')
plotting.plot_surf_roi(fsaverage['infl_right'], roi_map=col_fsa[int(col_fsa.shape[0]/2):],
hemi='right', view='lateral', bg_map=fsaverage['sulc_right'],
bg_on_data=True, darkness=0.5, cmap=ic_cmap, colorbar = False)
legend_elements = [
Patch(facecolor=color_shades3[0], label='Top 5%'),
# Patch(facecolor=color_shades3[1], label='Top 10%'),
Patch(facecolor=color_shades3[1], label='Top 15%'),
# Patch(facecolor=color_shades3[3], label='Top 20%'),
Patch(facecolor=color_shades3[2], label='Top 25%')
]
plt.legend(handles=legend_elements, loc='upper center', bbox_to_anchor=(0.5, 1.05),
title='Percentile Ranges', ncol=len(legend_elements), frameon=False)
plt.show()
data_ic = data_ic.reset_index(drop=True)
data_degree = data_degree.reset_index(drop=True)
#df_combined = pd.concat([data_ic, data_cmrglc, data_degree], axis=1)
df_combined = data_ic
#df_combined = data_cmrglc
#df_combined = data_degree.iloc[:,0:3]
df_combined.replace('Posterior_Cingulate', 'PCC', inplace=True)
df_combined.replace('Somatosensory_and_Motor', 'SomSen_Mot', inplace=True)
df_combined.replace('Dorsolateral_Prefrontal', 'DLPFC', inplace=True)
df_combined.replace('Dorsal_Stream_Visual', 'Dors_Str_Vis', inplace=True)
def add_padding(text, padding=1):
return ' ' * padding + text + ' ' * padding
k = 1
i = 0
col1 = ['#00915C','#e4f1e9']
col2 = ['#E6AC00','#fff4e3']
col3 = ['#0199DD','#D6F3FF']
number_of_columns = len(df_combined.columns)
col_width = 0.5 / number_of_columns # Equal width for all columns, adjust as needed
col_widths = [col_width] * number_of_columns # List of column widths
#--------------------------------------------------------------------------------------------
#---------------------- regions with strongest MwC values -------------------------------
#--------------------------------------------------------------------------------------------
df_combined = data_ic
df_combined.replace('Posterior_Cingulate', 'PCC', inplace=True)
df_combined.replace('Somatosensory_and_Motor', 'SomSen_Mot', inplace=True)
df_combined.replace('Dorsolateral_Prefrontal', 'DLPFC', inplace=True)
df_combined.replace('Dorsal_Stream_Visual', 'Dors_Str_Vis', inplace=True)
col = col1
header_colors = [col[i], col[i],col[i],col[i]]#, col2[i], col2[i],col2[i], col3[i], col3[i], col3[i]]
row_colors = [col[k], col[k], col[k], col[k]]#, col2[k], col2[k], col2[k], col3[k], col3[k], col3[k]]
fig, ax = plt.subplots(figsize=(7, 8), dpi=300)
ax.axis('tight')
ax.axis('off')
the_table = ax.table(cellText=df_combined.values,
colLabels=df_combined.columns,
cellLoc='center',
loc='center',
colWidths=col_widths)
the_table.auto_set_font_size(False)
the_table.set_fontsize(8)
for col, header_color in enumerate(header_colors):
header_cell = the_table[(0, col)]
header_cell.set_facecolor(header_color)
header_cell.set_text_props(color='white', weight='bold', size=8)
header_cell.set_height(0.3)
for col, row_color in enumerate(row_colors):
for row in range(1, len(df_combined) + 1):
body_cell = the_table[(row, col)]
body_cell.set_facecolor(row_color)
body_cell.set_text_props(color='black', weight='bold', size=6)
body_cell.set_height(0.1)
plt.savefig(f'../results/Figures/IC_hubs_table.png', dpi=300, bbox_inches='tight')
plt.tight_layout()
plt.show()
#--------------------------------------------------------------------------------------------
#---------------------- regions with strongest Degree values -------------------------------
#--------------------------------------------------------------------------------------------
df_combined = data_degree
df_combined.replace('Posterior_Cingulate', 'PCC', inplace=True)
df_combined.replace('Somatosensory_and_Motor', 'SomSen_Mot', inplace=True)
df_combined.replace('Dorsolateral_Prefrontal', 'DLPFC', inplace=True)
df_combined.replace('Dorsal_Stream_Visual', 'Dors_Str_Vis', inplace=True)
df_combined.replace('Paracentral_Lobular_and_Mid_Cingulate', 'PCL and MCC', inplace=True)
df_combined.replace('Inferior_', 'PCL and MCC', inplace=True)
col = col3
header_colors = [col[i], col[i],col[i],col[i]]#, col2[i], col2[i],col2[i], col3[i], col3[i], col3[i]]
row_colors = [col[k], col[k], col[k], col[k]]#, col2[k], col2[k], col2[k], col3[k], col3[k], col3[k]]
fig, ax = plt.subplots(figsize=(7,8), dpi=300)
ax.axis('tight')
ax.axis('off')
# Create the table with specified column widths
the_table = ax.table(cellText=df_combined.values,
colLabels=df_combined.columns,
cellLoc='center',
loc='center',
colWidths=col_widths)
the_table.auto_set_font_size(False)
the_table.set_fontsize(8) # Set the font size for all cells
# Set the background color and text properties for header cells
for col, header_color in enumerate(header_colors):
header_cell = the_table[(0, col)]
header_cell.set_facecolor(header_color)
header_cell.set_text_props(color='white', weight='bold', size=8)
header_cell.set_height(0.3) # Visually increase the header cell height
# Set the background color and text properties for the rest of the cells
for col, row_color in enumerate(row_colors):
for row in range(1, len(df_combined) + 1): # Start from 1 to skip the header
body_cell = the_table[(row, col)]
body_cell.set_facecolor(row_color)
body_cell.set_text_props(color='black', weight='bold', size=6)
body_cell.set_height(0.1) # Set a consistent height for the rest of the cells
plt.savefig(f'../results/Figures/DC_hubs_table.png', dpi=300, bbox_inches='tight')
plt.tight_layout()
plt.show()
# %% [markdown]
# ## Fig S2: external data (MwC)
# %%
COLOR1 = "#859F81"
COLOR2 = '#3C9A62'
session1 = "m1"
session2 = "m2"
session3 = "ZU"
LIMB = 'without'
data_path = f'../data/processed/'
with open(os.path.join(data_path,f'matrices_data_{session1}_{LIMB}_limbic.pkl'), 'rb') as file:
data_input_ex1 = pickle.load(file)
with open(os.path.join(data_path,f'matrices_data_{session2}_{LIMB}_limbic.pkl'), 'rb') as file:
data_input_ex2 = pickle.load(file)
with open(os.path.join(data_path,f'matrices_data_{session3}_{LIMB}_limbic.pkl'), 'rb') as file:
data_input_ex3 = pickle.load(file)
edata_medianallsub_ex1 = data_input_ex1['edata_medianallsub']
edata_medianallsub_ex2 = data_input_ex2['edata_medianallsub']
edata_medianallsub_ex3 = data_input_ex3['edata_medianallsub']
sub_size1 = edata_medianallsub_ex1.shape[1]
sub_size2 = edata_medianallsub_ex2.shape[1]
sub_size3 = edata_medianallsub_ex3.shape[1]
DATA_avg_ex1 = pd.read_csv(os.path.join(data_path,f'DATA_avg_woSC_{session1}_{LIMB}_limb.csv'))
with open(os.path.join(data_path,f'data_single_sub_woSC_{session1}_{LIMB}_limb.pkl'), 'rb') as file:
data_single_sub_ex1 = pickle.load(file)
DATA_avg_ex2 = pd.read_csv(os.path.join(data_path,f'DATA_avg_woSC_{session2}_{LIMB}_limb.csv'))
with open(os.path.join(data_path,f'data_single_sub_woSC_{session2}_{LIMB}_limb.pkl'), 'rb') as file:
data_single_sub_ex2 = pickle.load(file)
DATA_avg_ex3 = pd.read_csv(os.path.join(data_path,f'DATA_avg_wSC_{session3}_{LIMB}_limb.csv'))
with open(os.path.join(data_path,f'data_single_sub_wSC_{session3}_{LIMB}_limb.pkl'), 'rb') as file:
data_single_sub_ex3 = pickle.load(file)
ICallsub_w_ex1 = zscore(data_single_sub_ex1['ICallsub_w'], nan_policy='omit')
ICallsub_w_ex2 = zscore(data_single_sub_ex2['ICallsub_w'], nan_policy='omit')
ICallsub_w_ex3 = zscore(data_single_sub_ex3['ICallsub_w'], nan_policy='omit')
#--------------------------------------------------------------------------------------------
#------------------ SPatial Autocorrelation -----------------------
#--------------------------------------------------------------------------------------------
LIMB = 'without'
niter = 10
limb_ind = np.array([88, 90, 92, 93, 110, 118, 120, 122, 131, 135, 165, 166, 172, 268, 270, 272, 273, 290, 298, 300, 302, 311, 314, 315, 346, 352])-1
mode = "IC"
test_stat_ex1 ,surrogate_brainmap_corrs_ex1,sa_corrected_p_value_ex1,spatially_naive_p_value_ex1 = Spatial_AC(DATA_avg_ex1, mode, LIMB, niter )
test_stat_ex2 ,surrogate_brainmap_corrs_ex2,sa_corrected_p_value_ex2,spatially_naive_p_value_ex2 = Spatial_AC(DATA_avg_ex2, mode, LIMB, niter )
test_stat_ex3 ,surrogate_brainmap_corrs_ex3,sa_corrected_p_value_ex3,spatially_naive_p_value_ex3 = Spatial_AC(DATA_avg_ex3, mode, LIMB, niter )
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
sns.kdeplot(
surrogate_brainmap_corrs_ex1,
color=COLOR1,
fill=False,
linewidth=2.5,
alpha=0.7,
label="Data Ex1",
ax=ax
)
ax.axvline(
test_stat_ex1,
ymin=0, ymax=0.96,
color=COLOR1,
linestyle='dashed',
lw=3
)
sns.kdeplot(
surrogate_brainmap_corrs_ex2,
color=COLOR2,
fill=False,
linewidth=2.5,
alpha=0.7,
label="Data Ex2",
ax=ax
)
ax.axvline(
test_stat_ex2,
ymin=0, ymax=0.96,
color=COLOR2,
linestyle='dashed',
lw=3
)
sns.kdeplot(
surrogate_brainmap_corrs_ex3,
color=COLOR,
fill=False,
linewidth=2.5,
alpha=0.7,
label="Data Ex3",
ax=ax
)
ax.axvline(
test_stat_ex3,
ymin=0, ymax=0.96,
color=COLOR,
linestyle='dashed',
lw=3
)
ax.set_xticks(np.arange(-1, 1.1, 0.5))
ax.set_xlim(-0.8, 0.8)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.tick_params(axis='both', which='major', labelsize=FONT)
plt.xlabel('Correlation', fontsize=FONT)
plt.ylabel('Density', fontsize=FONT)
plt.show()
data_ex1_z = pd.DataFrame({'pet_avg':DATA_avg_ex1['pet_avg'], 'ICallsub_w_avg':zscore(DATA_avg_ex1['ICallsub_w_avg'], nan_policy='omit'),
'degallsub_w_avg':zscore(DATA_avg_ex1['degallsub_w_avg'], nan_policy='omit')})
data_ex2_z = pd.DataFrame({'pet_avg':DATA_avg_ex2['pet_avg'], 'ICallsub_w_avg':zscore(DATA_avg_ex2['ICallsub_w_avg'], nan_policy='omit'),
'degallsub_w_avg':zscore(DATA_avg_ex2['degallsub_w_avg'], nan_policy='omit')})
data_ex3_z = pd.DataFrame({'pet_avg':DATA_avg_ex3['pet_avg'], 'ICallsub_w_avg':zscore(DATA_avg_ex3['ICallsub_w_avg'], nan_policy='omit'),
'degallsub_w_avg':zscore(DATA_avg_ex3['degallsub_w_avg'], nan_policy='omit')})
corr_IC_w_E_ex1, p_IC_w_E_ex1 = pcor(data_ex1_z.ICallsub_w_avg, data_ex1_z.pet_avg)
corr_IC_w_deg_ex1, p_IC_w_deg_ex1 = pcor(data_ex1_z.ICallsub_w_avg, data_ex1_z.degallsub_w_avg)
corr_IC_w_E_ex2, p_IC_w_E_ex2 = pcor(data_ex2_z.ICallsub_w_avg, data_ex2_z.pet_avg)
corr_IC_w_deg_ex2, p_IC_w_deg_ex2 = pcor(data_ex2_z.ICallsub_w_avg, data_ex2_z.degallsub_w_avg)
corr_IC_w_E_ex3, p_IC_w_E_ex3 = pcor(data_ex3_z.ICallsub_w_avg, data_ex3_z.pet_avg)
corr_IC_w_deg_ex3, p_IC_w_deg_ex3 = pcor(data_ex3_z.ICallsub_w_avg, data_ex3_z.degallsub_w_avg)
corr_IC_E_allsub_ex1 = np.zeros((sub_size1))
p_IC_E_allsub_ex1 = np.zeros((sub_size1))
corr_IC_E_allsub_ex2 = np.zeros((sub_size2))
p_IC_E_allsub_ex2 = np.zeros((sub_size2))
corr_IC_E_allsub_ex3 = np.zeros((sub_size3))
p_IC_E_allsub_ex3 = np.zeros((sub_size3))
for j in range(sub_size1):
corr_IC_E_allsub_ex1[j], p_IC_E_allsub_ex1[j] = pcor(ICallsub_w_ex1[:, j], edata_medianallsub_ex1[:, j])
for j in range(sub_size2):
corr_IC_E_allsub_ex2[j], p_IC_E_allsub_ex2[j] = pcor(ICallsub_w_ex2[:, j], edata_medianallsub_ex2[:, j])
for j in range(sub_size3):
corr_IC_E_allsub_ex3[j], p_IC_E_allsub_ex3[j] = pcor(ICallsub_w_ex3[:, j], edata_medianallsub_ex3[:, j])
#----------------------------------------------------------------------------------------------------
# scatter IC_avg , CMRglc_avg / group analysis >>> external data Vienna with two sessions (m1 and m2)
#----------------------------------------------------------------------------------------------------
g = sns.jointplot(
data=data_ex1_z, x='pet_avg', y='ICallsub_w_avg', kind="reg",
scatter_kws={'s': 50, 'alpha': 0.7, 'facecolors': COLOR1, 'edgecolors': COLOR1},
line_kws={'color': COLOR1, 'lw': 2},
color=COLOR1, label="session1", marginal_ticks=False
)
g.ax_marg_x.clear()
g.ax_marg_y.clear()
g.ax_marg_x.tick_params(
axis='both',
which='both',
bottom=False,
labelbottom=False,
left=False,
labelleft=False
)
g.ax_marg_y.tick_params(
axis='both',
which='both',
left=False,
labelleft=False,
bottom=False,
labelbottom=False
)
sns.regplot(
data=data_ex2_z, x='pet_avg', y='ICallsub_w_avg',
scatter_kws={'s': 50, 'alpha': 0.7, 'facecolors': COLOR2, 'edgecolors': COLOR2},
line_kws={'color': COLOR2, 'lw': 2},
color=COLOR2, ax=g.ax_joint, label="session2"
)
sns.kdeplot(data=data_ex1_z, x='pet_avg', ax=g.ax_marg_x, color=COLOR1, lw=2, label="Session 1", fill=False)
sns.kdeplot(data=data_ex2_z, x='pet_avg', ax=g.ax_marg_x, color=COLOR2, lw=2, label="Session 2", fill=False)
sns.kdeplot(data=data_ex1_z, y='ICallsub_w_avg', ax=g.ax_marg_y, color=COLOR1, lw=2, label="Session 1", fill=False)
sns.kdeplot(data=data_ex2_z, y='ICallsub_w_avg', ax=g.ax_marg_y, color=COLOR2, lw=2, label="Session 2", fill=False)
g.ax_marg_x.set_ylabel("")
g.ax_marg_y.set_xlabel("")
g.ax_joint.set_xlabel('Target Energy (CMRglc)', fontsize=FONT)
g.ax_joint.set_ylabel('MwC (z-score)', fontsize=FONT)
g.ax_joint.set_ylim([-3,3])
p_val1 = sa_corrected_p_value_ex1;
p_val2 = p_IC_w_E_ex1
if p_IC_w_E_smash < 0.001:
title1 = "corr1 = {:.2f} , p-smash < {:.3f}".format(float(corr_IC_w_E_ex1 ), float(0.001))
else:
title1 = "corr1 = {:.2f} , p-smash < {:.3f}".format(float(corr_IC_w_E_ex1 ), float(p_val1.item()))
if p_IC_w_E_smash < 0.001:
title2 = "corr2 = {:.2f} , p-smash < {:.3f}".format(float(corr_IC_w_E_ex2 ), float(0.001))
else:
title2 = "corr2 = {:.2f} , p-smash < {:.3f}".format(float(corr_IC_w_E_ex2 ), float(p_val2.item()))
g.ax_joint.set_title(
title1 + "\n"+ title2,
fontsize=FONT, va='baseline', y=0.85
)
plt.grid(True)
g.ax_joint.tick_params(axis='both', which='major', labelsize=FONT)
spine_width = 1.6
for spine in g.ax_joint.spines.values():
spine.set_linewidth(spine_width)
g.ax_joint.legend(fontsize=FONT-2, loc="lower right")
plt.savefig('../results/Figures/ext_IC_CMRglc_wien.png', dpi=300, bbox_inches='tight')
plt.show()
#----------------------------------------------------------------------------------------------------
# scatter IC_avg , CMRglc_avg / group analysis >>> external data TUM: closed eyes(ZU):
#----------------------------------------------------------------------------------------------------
plt.figure(figsize=(7,7))
g = sns.jointplot(data=data_ex3_z, x='pet_avg', y='ICallsub_w_avg', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': COLOR, 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':IC_color, 'edgecolors':IC_color} )
g.ax_marg_x.clear()
g.ax_marg_y.clear()
g.ax_marg_x.tick_params(
axis='both',
which='both',
bottom=False,
labelbottom=False,
left=False,
labelleft=False
)
g.ax_marg_y.tick_params(
axis='both',
which='both',
left=False,
labelleft=False,
bottom=False,
labelbottom=False
)
sns.kdeplot(data=data_ex3_z, x='pet_avg', ax=g.ax_marg_x, color=COLOR, lw=2, fill=False)
sns.kdeplot(data=data_ex3_z, y='ICallsub_w_avg', ax=g.ax_marg_y, color=COLOR, lw=2 ,fill=False)
title = f"corr = {corr_IC_w_E_ex3:.2f}, p-smash< 0.001"
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'MCC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Target Energy (CMRglc)', fontsize = FONT)
plt.ylabel('MwC (z-score)', fontsize = FONT)
plt.ylim([-4,4])
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.savefig('../results/Figures/ext_IC_CMRglc_ZU.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
# Box plot individual correlation between IC and E >>> external data Vienna with two sessions (m1 and m2)
#--------------------------------------------------------------------------------------------
group_colors = [COLOR1, COLOR2]
data = pd.DataFrame({
'corr': np.concatenate([corr_IC_E_allsub_ex1, corr_IC_E_allsub_ex2]),
'session': ['session1'] * sub_size1+ ['session2'] * sub_size2,
'pval' : np.concatenate([p_IC_E_allsub_ex1, p_IC_E_allsub_ex2])
})
fig, ax = plt.subplots(figsize=(3, 6))
plt.grid(True)
violin = sns.violinplot(data=data, x = 'session', y='corr', ax=ax, color = COLOR)
violin.collections[0].set_facecolor('none')
violin.collections[1].set_facecolor('none')
sns.stripplot(data=data[data['pval'] < 0.05], x='session', y='corr', hue='session', palette=group_colors, size=8, ax=ax, jitter=True, label="p < 0.05")
sns.stripplot(data=data[data['pval'] >= 0.05], x='session', y='corr', color='red', size=8, ax=ax, jitter=True, label="p >= 0.05")
ax.set_ylabel('Correlation MwC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.15,0.72])
plt.tick_params(axis='both', which='major', labelsize= FONT, rotation = 45)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig('../results/Figures/ext_box_IC_CMRglc_wien.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
# Box plot individual correlation between IC and E >>> external data TUM (closed eyes)
#--------------------------------------------------------------------------------------------
data = pd.DataFrame({'corr': corr_IC_E_allsub_ex3, 'pval': p_IC_E_allsub_ex3})
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
violin = sns.violinplot(data=data, y='corr', ax=ax, color = COLOR)
violin.collections[0].set_facecolor('none')
sns.stripplot(data=data[data['pval'] < 0.05], y='corr', color=IC_color , size = 8, ax=ax)
sns.stripplot(data=data[data['pval'] >= 0.05], y='corr', color='red' , size = 8, ax=ax)
ax.set_ylabel('Correlation MwC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.15,0.72])
plt.tick_params(axis='both', which='major', labelsize= FONT)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig('../results/Figures/ext_box_IC_CMRglc_zu.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S3: participation coefficient
# %%
'''
participation coefficent can be calculated with one of these methods:
'infomap', 'louvain', 'walktrap', 'label_prop'
'''
PC_COLOR = '#7DA3DA'
session = "AUF"
LIMB = "without"
thr_fc = 0.1
thr_sc = 0.1
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_with_limb.csv'
DATA_avg = pd.read_csv(data_avg_path)
file_path = f'../data/processed/matrices_data_{session}_with_limbic.pkl'
with open(file_path, 'rb') as file:
data_input = pickle.load(file)
edata_medianallsub = data_input['edata_medianallsub']
conn_mat = data_input['Pearconn_allsub']
SCconn_allsub = data_input['SCconn_allsub']
sub_size = SCconn_allsub.shape[2]
all_pc = ParticipationCoeff(thr_fc, thr_sc, conn_mat, SCconn_allsub, sub_size)
pc_allmethods_avg = {method: np.mean(pc_matrix, axis=1) for method, pc_matrix in all_pc.items()}
data_pc_avg = pd.DataFrame({'pet_avg': DATA_avg.pet_avg, 'PCallsub_w_avg': pc_allmethods_avg[method], 'ICallsub_w_avg': DATA_avg.ICallsub_w_avg})
ind = data_pc_avg['PCallsub_w_avg']>0
#data_pc_avg = pd.DataFrame({'pet_avg': DATA_avg.pet_avg, 'PCallsub_w_avg': bg.pc, 'ICallsub_w_avg': DATA_avg.ICallsub_w_avg})
corr_PC_w_E, p_PC_w_E = pcor(data_pc_avg.PCallsub_w_avg[ind], data_pc_avg.pet_avg[ind])
#--------------------------------------------------------------------------------------------
#--------------------------- spatial autocorrelation correction -----------------------
#--------------------------------------------------------------------------------------------
N_ROT = 1000
RNG_SEED = 42
session = 'AUF'
LIMB = 'with'
ind = data_pc_avg['PCallsub_w_avg']>0
map1 = data_pc_avg.PCallsub_w_avg
map2 = data_pc_avg.pet_avg
cmrglc =map1.to_numpy(dtype=float).ravel()
ic = map2.to_numpy(dtype=float).ravel()
#rho_ic, p_spin_ic, null_ic = run_spin_test(ic, cmrglc, N_ROT,RNG_SEED)
# print(p_spin_ic)
FONT = 22
COLOR = IC_color
null_r = np.asarray(null_ic, dtype=float)
null_r = null_r[np.isfinite(null_r)]
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
sns.kdeplot(
null_r,
color="#E0E0E0",
ax=ax,
fill=True,
linewidth=2.5
)
ax.axvline(rho_ic, 0, 0.9, color='#9CB9E8', linestyle='dashed', lw=3)
ax.set_xticks(np.arange(-1, 1.1, 0.5))
ax.set_ylim(0, 5)
ax.set_xlim(-0.8, 0.8)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.tick_params(axis='both', which='major', labelsize=FONT)
plt.xlabel('Correlation', fontsize=FONT)
plt.ylabel('Density', fontsize=FONT)
plt.show()
lower_bound = np.percentile(null_r, 5)
upper_bound = np.percentile(null_r, 95)
print(f"90% CI of null distribution: [{lower_bound:.2f}, {upper_bound:.2f}]")
#--------------------------------------------------------------------------------------------
#--------------------------- group level: scatter plot CMRglc and PC -----------------------------------
#--------------------------------------------------------------------------------------------
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data_pc_avg, x='pet_avg', y='PCallsub_w_avg', kind="reg", color=PC_COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':PC_COLOR, 'edgecolors':PC_COLOR} )
p_PC_w_E_smash = p_spin_ic
if p_PC_w_E_smash < 0.001:
title = "corr = {:.2f} , psmash< {:.0e}".format(float(corr_PC_w_E), float(0.001))
else:
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_PC_w_E), float(p_PC_w_E_smash))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'PC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Target Energy (CMRglc)', fontsize = FONT)
plt.ylabel('PC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- subject level: box plot CMRglc and PC -----------------------------------
#--------------------------------------------------------------------------------------------
corr_PC_E_allsub = np.zeros((sub_size))
p_PC_E_allsub = np.zeros((sub_size))
for j in range(sub_size):
corr_PC_E_allsub[j], p_PC_E_allsub[j] = pcor(all_pc[method][:, j], edata_medianallsub[:, j])
data = pd.DataFrame({'corr': corr_PC_E_allsub, 'pval': p_PC_E_allsub})
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
violin = sns.violinplot(data=data, y='corr', ax=ax, color = PC_COLOR)
violin.collections[0].set_facecolor('none')
sns.stripplot(data=data[data['pval'] < 0.05], y='corr', color=PC_COLOR , size = 8, ax=ax)
sns.stripplot(data=data[data['pval'] >= 0.05], y='corr', color='red' , size = 8, ax=ax)
ax.set_ylabel('Correlation PC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.15,0.72])
plt.tick_params(axis='both', which='major', labelsize= FONT)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- PC map on the brain surface -----------------------------------
#--------------------------------------------------------------------------------------------
#plot_node_surf(data_pc_avg, Pearconn_allsub, 'PCallsub_w_avg', CMAP_DC,0.35, [])
PC = np.array(data_pc_avg.PCallsub_w_avg).reshape(-1, 1)
PC_allrois = np.full((360, 1), np.nan)
ind = np.arange(360)
mask = np.ones(ind.shape , dtype = bool)
#mask[limb_ind] = False
ind = ind[mask]
PC_allrois[ind] = np.abs(PC)
fig = plt.figure(figsize=(5, 5),dpi=300)
fsaverage = datasets.fetch_surf_fsaverage()
V_min = 0
V_max = PC.max()
#V_max = 34
ticks = [V_min, 0.1, V_max]
col_fsa = parcel_to_surface(PC_allrois,'glasser_360_fsa5')
plotting.plot_surf_stat_map( fsaverage['infl_left'], stat_map = col_fsa[:int(col_fsa.shape[0]/2)],
title='PC, left hemisphere',
hemi='left', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_left'], bg_on_data=False,
cmap=CMAP_DC,colorbar=False)
plt.show()
col_fsa = parcel_to_surface(PC_allrois,'glasser_360_fsa5')
plotting.plot_surf_stat_map(fsaverage['infl_left'], stat_map = col_fsa[:int(col_fsa.shape[0]/2)],
title='PC, left hemisphere',
hemi='left', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_left'], bg_on_data = False,
cmap=CMAP_DC,colorbar=False)
plt.show()
col_fsa = parcel_to_surface(PC_allrois,'glasser_360_fsa5')
plotting.plot_surf_stat_map(fsaverage['infl_right'], stat_map = col_fsa[int(col_fsa.shape[0]/2):],
title='PC, right hemisphere',
hemi='right', vmin= V_min, vmax = V_max, view='lateral',
bg_map=fsaverage['sulc_right'], bg_on_data=False,
cmap=CMAP_DC,colorbar=False)
plt.show()
col_fsa = parcel_to_surface(PC_allrois,'glasser_360_fsa5')
plotting.plot_surf_stat_map(fsaverage['infl_right'], stat_map = col_fsa[int(col_fsa.shape[0]/2):],
title='PC, right hemisphere',
hemi='right', vmin= V_min, vmax = V_max, view='medial',
bg_map=fsaverage['sulc_right'], bg_on_data = False,
cmap=CMAP_DC,colorbar=True)
plt.show()
norm = mpl.colors.Normalize(vmin=V_min, vmax=V_max)
sm = mpl.cm.ScalarMappable(cmap=CMAP_DC, norm=norm)
sm.set_array([])
fig, ax = plt.subplots(figsize=(6, 1))
cbar = plt.colorbar(
sm,
cax = ax,
orientation = 'horizontal',
ticks = ticks,
format = "%.2f"
)
cbar.set_label("Participation Coefficient", fontsize=10)
plt.tight_layout()
plt.show()
# %% [markdown]
# ## Fig S4: MwC vs CMRglc for combined datasets
# %%
#--------------------------------------------------------------------------------------------
#------------------ read all sessions data -------------------------------
#-------------------------------------------------------------------------------------------
session = 'AUF'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_wSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_auf = pickle.load(file)
DATA_avg_auf = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_auf = pickle.load(file)
DATA_avg_pr_auf = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_auf = pickle.load(file)
session = 'ZU'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_wSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_zu = pickle.load(file)
DATA_avg_zu = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_zu = pickle.load(file)
DATA_avg_pr_zu = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_zu = pickle.load(file)
session = 'm1'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_woSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_woSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_woSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_woSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_m1 = pickle.load(file)
DATA_avg_m1 = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_m1 = pickle.load(file)
DATA_avg_pr_m1 = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_m1 = pickle.load(file)
session = 'm2'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_woSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_woSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_woSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_woSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_m2 = pickle.load(file)
DATA_avg_m2 = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_m2 = pickle.load(file)
DATA_avg_pr_m2 = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_m2 = pickle.load(file)
#--------------------------------------------------------------------------------------------
#------------------ concatenated all sessions data -------------------------------
#-------------------------------------------------------------------------------------------
df1 = zscore(data_single_sub_auf['ICallsub_w'], nan_policy='omit')
df2 = zscore(data_single_sub_zu['ICallsub_w'], nan_policy='omit')
df3 = zscore(data_single_sub_m1['ICallsub_w'], nan_policy='omit')
df4 = zscore(data_single_sub_m2['ICallsub_w'], nan_policy='omit')
ICallsub_allsessions = np.hstack((df1, df2, df3, df4))
#ICallsub_allsessions = df3
ICallsub_allsessions = pd.DataFrame(ICallsub_allsessions)
df1 = zscore(data_single_sub_pr_auf['degallsub_w'], nan_policy='omit')
df2 = zscore(data_single_sub_pr_zu['degallsub_w'], nan_policy='omit')
df3 = zscore(data_single_sub_pr_m1['degallsub_w'], nan_policy='omit')
df4 = zscore(data_single_sub_pr_m2['degallsub_w'], nan_policy='omit')
degallsub_allsessions = np.hstack((df1, df2, df3, df4))
#degallsub_allsessions = df3
degallsub_allsessions = pd.DataFrame(degallsub_allsessions)
df1 = zscore(data_input_auf['edata_medianallsub'], nan_policy='omit')
df2 = zscore(data_input_zu['edata_medianallsub'], nan_policy='omit')
df3 = zscore(data_input_m1['edata_medianallsub'], nan_policy='omit')
df4 = zscore(data_input_m2['edata_medianallsub'], nan_policy='omit')
edata_allsub_allsessions = np.hstack((df1, df2, df3, df4))
#edata_allsub_allsessions = df3
edata_allsub_allsessions = pd.DataFrame(edata_allsub_allsessions)
ICallsub_allsessions_avg = np.nanmean(ICallsub_allsessions, axis = 1)
degallsub_allsessions_avg = np.nanmean(degallsub_allsessions, axis = 1)
edata_allsub_allsessions_avg = np.nanmean(edata_allsub_allsessions, axis = 1)
#--------------------------------------------------------------------------------------------
#------------------ concatenated all sessions data -------------------------------
#-------------------------------------------------------------------------------------------
IC = ICallsub_allsessions_avg
DC = degallsub_allsessions_avg
CMRglc = edata_allsub_allsessions_avg
IC = np.array(IC).flatten()
DC = np.array(DC).flatten()
CMRglc = np.array(CMRglc).flatten()
alldata_avg = pd.DataFrame({'IC':IC, 'DC':DC, 'pet':CMRglc})
r_IC_CMRglc,p_IC_CMRglc = pcor(IC,CMRglc)
print(r_IC_CMRglc )
r_DC_CMRglc , p_DC_CMRglc= pcor(DC,CMRglc)
print(r_DC_CMRglc)
r_IC_DC , p_IC_DC = pcor(IC,DC)
print(r_IC_DC)
#--------------------------------------------------------------------------------------------
#------------------ scatter IC_avg , CMRglc_avg / group analysis ----------------------------
#--------------------------------------------------------------------------------------------
COLOR = IC_color
plt.figure(figsize=(6,4))
g = sns.jointplot(data=alldata_avg, x='pet', y='IC', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "corr = {:.2f} , p-val = {:.0e}".format(float(r_IC_CMRglc), float(p_IC_CMRglc,))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'MCC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Target Energy (CMRglc)', fontsize = FONT)
plt.ylabel('MwC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.savefig('../results/Figures/all_scatter_IC_CMRglc_.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------ Box plot individual correlation between IC and E -----------------------
#--------------------------------------------------------------------------------------------
sub_size = np.size(edata_allsub_allsessions, axis =1)
corr_IC_E_allsub = np.zeros((sub_size))
p_IC_E_allsub = np.zeros((sub_size))
for j in range(sub_size):
corr_IC_E_allsub[j], p_IC_E_allsub[j] = pcor(edata_allsub_allsessions.iloc[:,j], ICallsub_allsessions.iloc[:,j])
data = pd.DataFrame({'corr': corr_IC_E_allsub, 'pval': p_IC_E_allsub})
fig, ax = plt.subplots(figsize=(3, 7))
plt.grid(True)
violin = sns.violinplot(data=data, y='corr', ax=ax, color = COLOR)
violin.collections[0].set_facecolor('none')
sns.stripplot(data=data[data['pval'] < 0.05], y='corr', color=COLOR , size = 8, ax=ax)
sns.stripplot(data=data[data['pval'] >= 0.05], y='corr', color='red' , size = 8, ax=ax)
ax.set_ylabel('Correlation MwC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.15,0.72])
plt.tick_params(axis='both', which='major', labelsize= FONT)
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig('../results/Figures/all_box_IC_CMRglc_.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S5-a: control analysis
# %%
#--------------------------------------------------------------------------------------------
#---- same weights for one target (keep the number of connection same as the real data)------
#--------------------------------------------------------------------------------------------
SCmat = SCconn_allsub
thr_sc = 0.2
SCmat_th = np.zeros((SCmat.shape[0], SCmat.shape[1], sub_size))
for i in range(sub_size):
SCmat_th[:, :, i] = threshold_proportional(SCmat[:, :, i], thr_sc)
SC_mask = SCmat_th.copy()
SC_mask[SC_mask > 0] = 1
Conn_Mat = SC_mask
# data_avg_sur2, data_single_sub_sur2 = IC_calculation(Conn_Mat, edata_medianallsub_rem, SCconn_allsub, thr_sc, thr_FC, sub_size, scmask, nrois_rem,net_label)
# DATA_avg_sur2 = pd.concat([data_avg_sur2 , net_label], axis =1)
corr_IC_E_allsub = np.zeros((sub_size))
p_IC_E_allsub = np.zeros((sub_size))
for j in range(sub_size):
corr_IC_E_allsub[j], p_IC_E_allsub[j] = pcor(data_single_sub_sur2["ICallsub_w"][:, j], edata_medianallsub_rem[:, j])
organg = [255/255,164/255,0/255]
corr_IC_w_E, p_IC_w_E = pcor(DATA_avg_sur2.ICallsub_w_avg, DATA_avg_sur2.pet_avg)
corr_IC_w_deg, p_IC_w_deg = pcor(DATA_avg_sur2.ICallsub_w_avg, DATA_avg_sur2.degallsub_w_avg)
g = sns.jointplot(data = DATA_avg_sur2 , x = 'pet_avg' , y ='ICallsub_w_avg' , kind="reg", color="#C83D64",
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':"#E6E6E6", 'edgecolors':"#787B76"} )
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_w_E) , float(p_IC_w_E ))
g.fig.suptitle(f"corr = {corr_IC_w_E:.2f} , p-val = {p_IC_w_E :.0e}", fontsize=FONT, va='baseline', y=0.8)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.xlabel('Energy of target', fontsize = FONT)
plt.ylabel('MwC (CMRglc-only)', fontsize = FONT)
plt.grid(True)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.savefig('../results/Figures/control_scatter_CMRglc_only.png', dpi=300, bbox_inches='tight')
plt.show()
#------------------ Box plot individual correlation between IC and E -----------------------
print(p_IC_E_allsub)
data = pd.DataFrame({'corr': corr_IC_E_allsub, 'pval': p_IC_E_allsub})
fig, ax = plt.subplots(figsize=(3, 6))
plt.grid(True)
violin = sns.violinplot(data=data, y='corr', ax=ax, color = '#FBE5CF')
sns.stripplot(data=data[data['pval'] < 0.05], y='corr', color='#787B76' , size = 8, ax=ax)
sns.stripplot(data=data[data['pval'] >= 0.05], y='corr', color='red' , size = 8, ax=ax)
violin.collections[0].set_facecolor('none')
ax.set_ylabel('Correlation MwC and target\'s energy', fontsize = FONT)
ax.set_ylim([-0.1,0.72])
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig('../results/Figures/control_box_CMRglc_only.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#--------------------------- permutation test ----------------------------------------------
#--------------------------------------------------------------------------------------------
#1. permutation on average data (IC_avg , E_avg)
observed_correlation, observed_pvalue = pcor(DATA_avg.pet_avg, DATA_avg.ICallsub_w_avg)
n_permutations = 500
permuted_correlations = np.zeros(n_permutations)
for perm_idx in range(n_permutations):
permuted_e_avg_values = np.copy(DATA_avg.pet_avg)
np.random.shuffle(permuted_e_avg_values)
permuted_correlations[perm_idx], _ = pcor(DATA_avg.ICallsub_w_avg, permuted_e_avg_values)
#p_value_permutation = np.mean(np.abs(permuted_correlations) >= np.abs(observed_correlation))
corr_perm_avg = np.mean(permuted_correlations)
# Print the results
print("Permutation Test Results:")
print(f"Observed Correlation: {observed_correlation}")
print(f"Observed p-value: {observed_pvalue}")
print(f"Average perumuted correlation : {corr_perm_avg }")
sns.displot(permuted_correlations, kde=True, color='#787B76')
#sns.xlabel('Correlation IE and target\'s energy')
plt.savefig('../results/Figures/control_RandomMatch_CMRglc_only.png', dpi=300, bbox_inches='tight')
plt.show()
#2. permutation on individual data (IC , E)
permuted_correlations = np.zeros((n_permutations, sub_size))
for i in range(sub_size):
E_sub = edata_medianallsub_rem[:,i]
IC_sub = ICallsub_w[:,i]
for perm_idx in range(n_permutations):
permuted_E_values = np.copy(E_sub)
np.random.shuffle(permuted_E_values)
permuted_correlations[perm_idx,i], _ = pcor(IC_sub, permuted_E_values)
#sns.boxplot(permuted_correlations)
#plt.show()
#---------------------- Ridge plot of each subject histogram -------------------------------
num_columns = permuted_correlations.shape[1]
reshaped_array = np.vstack((permuted_correlations.flatten(), np.tile(np.arange(num_columns), permuted_correlations.shape[0]))).T
df_reshaped = pd.DataFrame(reshaped_array, columns=['Correlation Value', 'Subject'])
df_filtered = df_reshaped
sns.set_theme(style="white", rc={"axes.facecolor": (0, 0, 0, 0), 'axes.linewidth':1})
palette = "Spectral"
g = sns.FacetGrid(df_filtered, palette=palette, row="Subject", hue="Subject", aspect=9, height=0.3)
g.map_dataframe(sns.kdeplot, x='Correlation Value', fill=True, alpha=0.8)
g.map_dataframe(sns.kdeplot, x='Correlation Value', color='black', linewidth= 0.7)
def label(x, color, label):
ax = plt.gca()
ax.text(0, .2, label, color='black', fontsize=10,
ha="left", va="center", transform=ax.transAxes)
g.fig.subplots_adjust(hspace=-.5)
g.set_titles("")
g.set(yticks=[], ylabel=None)
g.despine(left=True)
for ax in g.axes.flat:
ax.axvline(0, color='gray', linestyle='--', linewidth=1)
g.despine(left=True)
plt.savefig('../results/Figures/control_RandomMatch_single_CMRglc_only.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S5-b: Sensitivity analysis
# %%
#--------------------------------------------------------------------------------------------
#------------------------- sensitivity to sc threshold -----------------------------------
#--------------------------------------------------------------------------------------------
plt.grid(True)
edata_medianallsub_rem = edata_medianallsub
nrois = 360
nrois_rem = nrois - len(limb_ind)
sub_size = MIconn_allsub.shape[2]
MIconn_allsub[MIconn_allsub < 0.01]=0
Conn_Mat = MIconn_allsub
thr_FC = 1
scmask = 1
thr_sc = np.linspace(0.01,1 , 20 )
corr_scthr = np.zeros((len(thr_sc)))
p_scthr = np.zeros((len(thr_sc)))
for i , thr in enumerate(thr_sc):
IC = only_IC_avg(Conn_Mat, edata_medianallsub_rem, SCconn_allsub, thr, thr_FC, sub_size, scmask, nrois_rem)
corr_scthr[i] , p_scthr[i] = pcor(IC, DATA_avg.pet_avg)
sc = plt.scatter(thr_sc, corr_scthr, s=100, c=-np.log10(p_scthr), cmap='plasma')
plt.colorbar(sc, label='-log10(p-value)')
plt.ylim([0.5,0.7])
plt.xlabel('SC threshold', fontsize = FONT)
plt.ylabel('Corr MwC and target\'s energy', fontsize = FONT)
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.savefig('../results/Figures/SC_threshold.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------------- sensitivity to fc threshold -----------------------------------
#--------------------------------------------------------------------------------------------
plt.grid(True)
edata_medianallsub_rem = edata_medianallsub
nrois = 360
nrois_rem = nrois - len(limb_ind)
sub_size = MIconn_allsub.shape[2]
MIconn_allsub[MIconn_allsub < 0.0001]=0
Conn_Mat = MIconn_allsub
thr_sc = 1
scmask = 0
thr_fc = np.linspace(0.1,1 , 20 )
corr_fcthr = np.zeros((len(thr_fc)))
p_fcthr = np.zeros((len(thr_fc)))
for i , thr in enumerate(thr_fc):
IC = only_IC_avg(Conn_Mat, edata_medianallsub_rem, SCconn_allsub, thr_sc, thr, sub_size, scmask, nrois_rem)
corr_fcthr[i] , p_fcthr[i] = pcor(IC, DATA_avg.pet_avg)
sc = plt.scatter(thr_fc, corr_fcthr, s=100, c=-np.log10(p_fcthr), cmap='plasma')
plt.colorbar(sc, label='-log10(p-value)')
plt.ylim([0.5,0.7])
plt.xlabel('MI threshold', fontsize = FONT)
plt.ylabel('Corr MwC and target\'s energy', fontsize = FONT)
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.savefig('../results/Figures/MI_threshold.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------- sensitivity to region size (number of voxels) -----------------------
#--------------------------------------------------------------------------------------------
COLOR = [227/255 , 108/255, 138/255]
num_voxels_avg = np.mean(num_voxels , axis = 1)
corr_IC_w_vox, p_IC_w_vox = pcor(DATA_avg.ICallsub_w_avg, num_voxels_avg)
plt.figure(figsize=(6,4))
g = sns.jointplot(x = num_voxels_avg , y = DATA_avg.ICallsub_w_avg, kind="reg", color="#547BC9",
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':"#FDC827", 'edgecolors':"#FDC827"} )
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_w_vox) , float(p_IC_w_vox))
g.fig.suptitle(f"corr = {corr_IC_w_vox:.2f} , p-val = {p_IC_w_vox:.0e}", fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Actual Target Energy', 'MCC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.xlabel('Target size (#voxels)', fontsize = FONT)
plt.ylabel('MwC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.grid(True)
plt.savefig('../results/Figures/Target_size.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#-------------------------- 3. sensitivity to source numbers --------------------------------
#--------------------------------------------------------------------------------------------
corr_IC_w_source, p_IC_w_source = pcor(DATA_avg.ICallsub_w_avg, DATA_avg.degallsub_b_avg)
#sns.lmplot(data = DATA_avg , x = 'pet_avg' , y ='ICallsu_w_avg' , hue = 'yeo_7_nw' )
sns.scatterplot(data = DATA_avg , x = 'degallsub_b_avg' , y ='ICallsub_w_avg' , color = 'Orange')
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_w_source) , float(p_IC_w_source))
plt.title(title , fontdict={'fontsize': 15})
plt.xlabel('Source Numbers')
plt.ylabel('MwC')
plt.show()
# %% [markdown]
# ## Fig S6: replication spatial correlation using Julich parcellation
# %%
#--------------------------------------------------------------------------------------------
#------------------ remove limbic and subcortical rois ----------------------------
#--------------------------------------------------------------------------------------------
atlas = '../data/external/JulichBrainAtlas_fsl_mni_3mm.nii.gz'
atlas_data = nib.load(atlas).get_fdata()
all_labels = np.unique(atlas_data)
all_labels = all_labels[all_labels > 0]
counts_remaining = np.array([np.sum(atlas_data == lab) for lab in all_labels])
small_rois = all_labels[counts_remaining <= 10] # exclude small regions
removed_rois = [17, 130, 120] # these three rois have more than 75% overlap with MMP limbic's rois
subcort_ids= [5, 9, 12, 15, 16,17, 18, 23, 24, 31, 38, 39, 44, 45, 47, 51, 52,
56, 61, 63, 66, 68, 74, 75, 79, 81, 87, 88, 92, 93, 103, 109,
113, 115, 119, 120, 125, 131, 136, 137, 138, 141, 146, 149, 151,
153, 154, 156, 166, 168, 169, 172, 173, 174, 176, 178, 181, 184,
185, 187, 188, 190, 191, 193, 194, 202, 212, 216, 219, 222, 223,224,
225, 230, 231, 238, 245, 246, 251, 252, 256, 259, 260, 263, 268,
270, 271, 275, 276, 280, 281, 291, 296, 300, 302, 306, 307, 312,
318, 320, 322, 323, 326, 331, 334, 336, 339, 340, 341, 343, 344,
353, 355, 356, 358, 361, 364, 365, 367, 371, 372, 373, 375, 378,
381, 382, 384,385, 388, 391, 392, 394, 395, 397, 398, 400, 401, 409]
rem_ind = np.union1d(removed_rois, subcort_ids)
rem_ind = np.union1d(rem_ind, small_rois)
roi_ids = np.unique(atlas_data[atlas_data != 0]).astype(np.int32)
roi_ids = np.setdiff1d(roi_ids, np.array(rem_ind, dtype=np.int32))
roi_ids = np.sort(roi_ids)
#print(roi_ids.shape)
#--------------------------------------------------------------------------------------------
#-------------- read functional connectivity and other inputs ----------------
#--------------------------------------------------------------------------------------------
session = "AUF"
LIMB = "without"
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic_Julich.pkl'
with open(file_path, 'rb') as file:
data_input_ju = pickle.load(file)
num_voxels_ju = data_input_ju['num_voxels']
MIconn_allsub_ju = data_input_ju['MIconn_allsub']
Pearconn_allsub_ju = data_input_ju['Pearconn_allsub']
edata_medianallsub_ju = data_input_ju['edata_medianallsub']
SCconn_allsub_ju = np.ones(Pearconn_allsub_ju.shape)
nrois_rem = MIconn_allsub_ju.shape[0]
sub_size = MIconn_allsub_ju.shape[2]
#--------------------------------------------------------------------------------------------
#-------------- calculate and save IC and other inputs ----------------
#--------------------------------------------------------------------------------------------
data_avg_pr_path = f'../data/processed/DATA_avg_woSC_pr_{session}_{LIMB}_limb_Julich.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_woSC_pr_{session}_{LIMB}_limb_Julich.pkl'
data_avg_path = f'../data/processed/DATA_avg_woSC_{session}_{LIMB}_limb_Julich.csv'
data_sub_path = f'../data/processed/data_single_sub_woSC_{session}_{LIMB}_limb_Julich.pkl'
MIconn_allsub_ju[MIconn_allsub_ju < 0.01] = 0
Conn_Mat = MIconn_allsub_ju
thr_FC = 1
scmask = 0
thr_sc = 0.1
[data_avg, data_single_sub] = IC_calculation(Conn_Mat, edata_medianallsub_ju, SCconn_allsub_ju, thr_sc, thr_FC, sub_size, scmask, nrois_rem)
DATA_avg = data_avg
DATA_avg.insert(0, "roi_ids_rem", roi_ids)
DATA_avg.to_csv(data_avg_path, index=False)
with open(data_sub_path, 'wb') as file:
pickle.dump(data_single_sub, file, protocol=pickle.HIGHEST_PROTOCOL)
# DATA_avg = pd.read_csv(data_avg_path)
# with open(data_sub_path, 'rb') as file:
# data_single_sub = pickle.load(file)
ICallsub_w = data_single_sub['ICallsub_w']
ICallsub_b = data_single_sub['ICallsub_b']
degallsub_w = data_single_sub['degallsub_w']
degallsub_b = data_single_sub['degallsub_b']
AvgMIallsub = data_single_sub['AvgMIallsub']
#--------------------------------------------------------------------------------------------
#------------------ scatter IC_avg , CMRglc_avg / group analysis ----------------------------
#--------------------------------------------------------------------------------------------
DATA_avg_new = DATA_avg.copy()
# DATA_avg_new = DATA_avg_new[DATA_avg_new ["pet_avg"]>20]
# DATA_avg_new = DATA_avg_new[DATA_avg_new ["ICallsub_w_avg"]>28.5]
corr_IC_w_E, p_IC_w_E = pcor(DATA_avg_new.ICallsub_w_avg, DATA_avg_new.pet_avg)
COLOR = "#8DD084"
plt.figure(figsize=(6,4))
g = sns.jointplot(data=DATA_avg_new, x='pet_avg', y='ICallsub_w_avg', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
p_IC_w_E_smash = p_IC_w_E
if p_IC_w_E_smash < 0.0001:
title = "r = {:.2f} , p-val < {:.3f}".format(float(corr_IC_w_E), float(0.001))
else:
title = "r = {:.2f} , p-val = {:.3f}".format(float(corr_IC_w_E), float(p_IC_w_E_smash))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'MwC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.ylim([28,31.5])
plt.xlabel('Target Energy (CMRglc)', fontsize = FONT)
plt.ylabel('MwC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
#plt.savefig(f'../results/Figures/CMRglc_IC_scatter_Julich.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S7: Effect of registraion quality on spatial correlation
# %%
sub = [3,7,12,14,17,20,23,25,26,28,29,30,31,32,33,35,36,37,38]
session = "AUF"
LIMB = "without"
sub_size = len(sub)
data_sub_path = f'../data/processed/data_single_sub_wSC_{session}_{LIMB}_limb.pkl'
data_sub_pr_path = f'../data/processed/data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
with open(file_path, 'rb') as file:
data_input = pickle.load(file)
edata_medianallsub = data_input['edata_medianallsub']
with open(data_sub_path, 'rb') as file:
data_single_sub = pickle.load(file)
ICallsub_w = data_single_sub['ICallsub_w']
degallsub_w = data_single_sub['degallsub_w']
print(degallsub_w.shape)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr = pickle.load(file)
ICallsub_w_pr = data_single_sub_pr['ICallsub_w']
degallsub_w_pr = data_single_sub_pr['degallsub_w']
corr_IC_E_allsub = np.zeros(sub_size)
p_IC_E_allsub = np.zeros(sub_size)
corr_deg_w_E_allsub = np.zeros(sub_size)
p_deg_w_E_allsub = np.zeros(sub_size)
for j in range(sub_size):
corr_IC_E_allsub[j], p_IC_E_allsub[j] = pcor(ICallsub_w[:, j], edata_medianallsub[:, j])
corr_deg_w_E_allsub[j], p_deg_w_E_allsub[j] = pcor(degallsub_w[:, j], edata_medianallsub[:, j])
######
mni2func_qc = pd.read_csv( "../data/processed/MNI2FUNC_AUF_qc_corratio_steps.csv")
mni2pet_qc = pd.read_csv( "../data/processed/MNI2PET_AUF_qc_corratio_steps.csv")
data = pd.DataFrame({"mni2func_qc":mni2func_qc["final"], "mni2pet_qc":mni2pet_qc["final"],
"corr_IC_E_allsub":corr_IC_E_allsub, "corr_deg_w_E_allsub":corr_deg_w_E_allsub})
corr_IC_qc , p_IC_qc = pcor(mni2func_qc["final"], corr_IC_E_allsub)
print(corr_IC_qc , p_IC_qc)
FRONT = 20
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='corr_IC_E_allsub', y='mni2func_qc', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "r = {:.2f} , p = {:.2f}".format(float(corr_IC_qc), float(p_IC_qc))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'Degree', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('subject-level corr(MwC,CMRglc)', fontsize = FONT)
plt.ylabel('corratio(MNI to functional)', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig(f'../results/Figures/qc_mni2func_corr_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
corr_IC_qc_mni2pet , p_IC_qc_mni2pet= pcor(mni2pet_qc["final"], corr_IC_E_allsub)
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='corr_IC_E_allsub', y='mni2pet_qc', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "r = {:.2f} , p = {:.2f}".format(float(corr_IC_qc_mni2pet), float(p_IC_qc_mni2pet))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'Degree', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('subject-level corr(MwC,CMRglc)', fontsize = FONT)
plt.ylabel('corratio(MNI to PET)', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
spine_color = 'gray'
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_color(spine_color)
plt.savefig(f'../results/Figures/qc_mni2pet_corr_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S8: relationship between functional/microstructural gradient and MwC
# %% [markdown]
# ### functional gradient (Margulies 2016)
# %%
LIMB = "without"
session = "AUF"
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
grad1_region_rem= pd.read_csv('../data/external/fcgrad1_region_rem.csv', header = None)
grad1_region_rem = grad1_region_rem.to_numpy(dtype=float).ravel()
DATA_avg = pd.read_csv(data_avg_path)
DATA_avg_pr = pd.read_csv(data_avg_pr_path)
IC =DATA_avg['ICallsub_w_avg']
deg = DATA_avg_pr['degallsub_w_avg']
pet = DATA_avg['pet_avg']
cor_IC_grad , p_IC_grad = pearsonr(IC,grad1_region_rem )
cor_deg_grad , p_deg_grad = pearsonr(deg,grad1_region_rem )
print(cor_IC_grad , p_IC_grad)
print(cor_deg_grad , p_deg_grad)
data = pd.DataFrame({"IC":IC, "deg":deg, "pet":pet, "grad1":grad1_region_rem})
#--------------------------------------------------------------------------------------------
#------------------ scatter IC_avg , fc gradient group analysis ----------------------------
#--------------------------------------------------------------------------------------------
COLOR = 'gray'
FONT = 18
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='grad1', y='IC', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "r = {:.2f} , p = {:.3f}".format(float(cor_IC_grad), float(p_IC_grad))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('First gradient', 'MwC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('First gradient', fontsize = FONT)
plt.ylabel('MwC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
plt.savefig(f'../results/Figures/grad1_IC_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------ scatter deg_avg , first gradient ----------------------------
#--------------------------------------------------------------------------------------------
COLOR = 'gray'
FONT = 18
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='grad1', y='deg', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "r = {:.2f} , p = {:.3f}".format(float(cor_deg_grad), float(p_deg_grad))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('First gradient', 'DC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('First gradient', fontsize = FONT)
plt.ylabel('wDC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
plt.savefig(f'../results/Figures/grad1_DC_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ### microstructural gradient
# %%
LIMB = "without"
session = "AUF"
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
grad1_region_rem = pd.read_csv('../data/external/micro_grad1_region_rem.csv', header = None)
grad1_region_rem = grad1_region_rem.to_numpy(dtype=float).ravel()
DATA_avg = pd.read_csv(data_avg_path)
DATA_avg_pr = pd.read_csv(data_avg_pr_path)
IC =DATA_avg['ICallsub_w_avg']
deg = DATA_avg_pr['degallsub_w_avg']
pet = DATA_avg['pet_avg']
cor_IC_grad , p_IC_grad = pearsonr(IC,grad1_region_rem )
cor_deg_grad , p_deg_grad = pearsonr(deg,grad1_region_rem )
print(cor_IC_grad , p_IC_grad)
print(cor_deg_grad , p_deg_grad)
data = pd.DataFrame({"IC":IC, "deg":deg, "pet":pet, "grad1":grad1_region_rem})
#--------------------------------------------------------------------------------------------
#------------------ scatter IC_avg , CMRglc_avg / group analysis ----------------------------
#--------------------------------------------------------------------------------------------
COLOR = 'gray'
FONT = 18
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='grad1', y='IC', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "r = {:.2f} , p = {:.3f}".format(float(cor_IC_grad), float(p_IC_grad))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('First gradient', 'MwC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Functional gradient', fontsize = FONT)
plt.ylabel('MwC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
plt.savefig(f'../results/Figures/micro_grad_IC_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
#--------------------------------------------------------------------------------------------
#------------------ scatter IC_avg , CMRglc_avg / group analysis ----------------------------
#--------------------------------------------------------------------------------------------
COLOR = 'gray'
FONT = 18
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='grad1', y='deg', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "r = {:.2f} , p = {:.3f}".format(float(cor_deg_grad), float(p_deg_grad))
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('First gradient', 'DC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('Microstructural gradient', fontsize = FONT)
plt.ylabel('DC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
spine.set_edgecolor("gray")
plt.savefig(f'../results/Figures/micro_grad_DC_scatter.png', dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## Fig S10: MwC-Mito association
# %%
mito_img_CI = nib.load('PATH to CI.nii.gz')
mito_img_CII = nib.load('PATH to CII.nii.gz')
mito_img_CIV = nib.load('PATH to CIV.nii.gz')
mito_img_TRC = nib.load('PATH to TRC.nii.gz')
mito_img_MRC = nib.load('PATH to MRC.nii.gz')
mito_img_MitoD = nib.load('PATH to MitoD.nii.gz')
mito_data_CI = mito_img_CI.get_fdata()
mito_data_CII = mito_img_CII.get_fdata()
mito_data_CIV = mito_img_CIV.get_fdata()
mito_data_TRC = mito_img_TRC.get_fdata()
mito_data_MRC = mito_img_MRC.get_fdata()
mito_data_MitoD = mito_img_MitoD.get_fdata()
mmp_in_mni_6 = f'../data/external/mmp_in_mni_6.nii.gz'
atlas_mito = nib.load(mmp_in_mni_6)
atlas_data_mito = atlas_mito.get_fdata()
region_labels_mito = np.unique(atlas_data_mito)
region_labels_mito = region_labels_mito[region_labels_mito > 0]
n_regions = len(region_labels_mito)
mito_values = {
'CI': np.zeros(n_regions),
'CII': np.zeros(n_regions),
'CIV': np.zeros(n_regions),
'TRC': np.zeros(n_regions),
'MRC': np.zeros(n_regions),
'MitoD': np.zeros(n_regions),
'syn': np.zeros(n_regions),
}
for i, label in enumerate(region_labels_mito):
region_mask = atlas_data_mito == label
region_vals_CI = mito_data_CI[region_mask]
region_vals_CII = mito_data_CII[region_mask]
region_vals_CIV = mito_data_CIV[region_mask]
region_vals_TRC = mito_data_TRC[region_mask]
region_vals_MRC = mito_data_MRC[region_mask]
region_vals_MitoD = mito_data_MitoD[region_mask]
# Compute mean of nonzero, non-NaN voxels
mito_values['CI'][i] = np.nanmedian(region_vals_CI[region_vals_CI != 0])
mito_values['CII'][i] = np.nanmedian(region_vals_CII[region_vals_CII != 0])
mito_values['CIV'][i] = np.nanmedian(region_vals_CIV[region_vals_CIV != 0])
mito_values['TRC'][i] = np.nanmedian(region_vals_TRC[region_vals_TRC != 0])
mito_values['MRC'][i] = np.nanmedian(region_vals_MRC[region_vals_MRC != 0])
mito_values['MitoD'][i] = np.nanmedian(region_vals_MitoD[region_vals_MitoD != 0])
# delete limbic regions, if needed:
mito_values['CI'] = np.delete(mito_values['CI'], limb_ind , axis = 0)
mito_values['CII'] = np.delete(mito_values['CII'], limb_ind , axis = 0)
mito_values['CIV'] = np.delete(mito_values['CIV'], limb_ind , axis = 0)
mito_values['TRC'] = np.delete(mito_values['TRC'], limb_ind , axis = 0)
mito_values['MRC'] = np.delete(mito_values['MRC'], limb_ind , axis = 0)
mito_values['MitoD'] = np.delete(mito_values['MitoD'], limb_ind , axis = 0)
data = pd.DataFrame({'IC':(DATA_avg.ICallsub_w_avg), 'deg':(DATA_avg_pr.degallsub_w_avg), 'pet':(DATA_avg.pet_avg)
,'mito_CI':(mito_values['CI']), 'mito_CII':(mito_values['CII']), 'mito_CIV':(mito_values['CIV']), 'mito_TRC':(mito_values['TRC'])
, 'mito_MRC':(mito_values['MRC']), 'mito_MitoD':(mito_values['MitoD'])})
corr_IC_CI, p_IC_CI =spearmanr(data["mito_CI"], data["IC"])
corr_IC_CII, p_IC_CII =spearmanr(data["mito_CII"], data["IC"])
corr_IC_CIV, p_IC_CIV = spearmanr(data["mito_CIV"], data["IC"])
corr_IC_TRC, p_IC_TRC = spearmanr(data["mito_TRC"], data["IC"])
corr_IC_MRC, p_IC_MRC = spearmanr(data["mito_MRC"], data["IC"])
corr_IC_MitoD, p_IC_MitoD = spearmanr(data["mito_MitoD"], data["IC"])
COLOR = IC_color
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='IC', y='mito_MitoD', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_MitoD), float(p_IC_MitoD))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
g.set_axis_labels('Target Energy', 'MCC', fontsize=FONT)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('MwC', fontsize = FONT)
plt.ylabel('Mitochondrial Density', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.show()
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='IC', y='mito_MRC', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_MRC), float(p_IC_MRC))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('MwC', fontsize = FONT)
plt.ylabel('MRC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.show()
plt.figure(figsize=(6,4))
g = sns.jointplot(data=data, x='IC', y='mito_TRC', kind="reg", color=COLOR,
marginal_kws=dict(bins=15, fill=True, color='gray'),
line_kws={'color': 'black', 'lw': 1},
scatter_kws={'s': 85, 'alpha': 0.7, 'facecolors':COLOR, 'edgecolors':COLOR} )
title = "corr = {:.2f} , p-val = {:.0e}".format(float(corr_IC_TRC), float(p_IC_TRC))
#g.fig.suptitle(f"corr = {corr_IC_w_E:.2f}"", $Pval_{smash}$ < 0.001", fontsize=FONT, va='baseline', y=0.8)
g.fig.suptitle(title, fontsize=FONT, va='baseline', y=0.8)
plt.tick_params(axis='both', which='major', labelsize= FONT)
plt.grid(True)
plt.xlabel('MwC', fontsize = FONT)
plt.ylabel('TRC', fontsize = FONT)
ax = g.ax_joint
spine_width = 1.6
for spine in ax.spines.values():
spine.set_linewidth(spine_width)
plt.show()
# %% [markdown]
# ## Fig S11:CMRglc - networks
# %%
net_label_WithLimbic = pd.read_csv(f'../data/external/mmp2yeo7nw_mapping.csv')
net_label_WithLimbic.reset_index(drop=True, inplace=True)
network_to_number = {network: i + 1 for i, network in enumerate(net_label_WithLimbic['yeo_7_nw'].unique())}
net_label_WithLimbic['network_number'] = net_label_WithLimbic['yeo_7_nw'].map(network_to_number)
color = [(139, 19, 140),(1,131, 182),(51, 116, 32),
(222, 75, 82),(226, 55, 255),(239, 156, 60), (255, 254, 211)]
net_color_WithLimbic = [(r/255, g/255, b/255 , 1) for r,g,b in color]
network_palette = dict(zip(net_names, net_color_WithLimbic))
xtik = ['Vis','Som','Dors','Def', 'Sal' , 'Cont', 'Limbic']
net_names_WithLimbic = ['Vis', 'SomMot' ,'DorsAttn' ,'Default' ,'SalVentAttn', 'Cont', 'Limbic' ]
net_num_WithLimbic = net_label_WithLimbic["network_number"]
new_net_num_WithLimbic = np.tile(net_num.transpose(), (nrois, 1))
num_net = len(net_names_WithLimbic)
LIMB = "with"
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
DATA_avg_WithLimbic = pd.read_csv(data_avg_path)
edata = DATA_avg_WithLimbic.pet_avg
flat_values = edata.values
flat_networks = net_label_WithLimbic['yeo_7_nw'].values
df_box = pd.DataFrame({
'Value': flat_values,
'Network': flat_networks
})
plt.figure(figsize=(10, 6))
sns.set_style("whitegrid")
sns.boxplot(x='Network', y='Value', data=df_box, palette=network_palette)
plt.xticks(rotation=45)
plt.ylabel('Median CMRglc')
plt.xlabel('Network')
plt.title(f'Median CMRglc across Networks with Limbic')
plt.tight_layout()
plt.show()
# %% [markdown]
# # Statistical comparison DC, MCC
# %% [markdown]
# ## regression CMRglc = DC + MCC
# %%
#--------------------------------------------------------------------------------------------
#------------------ read all sessions data -------------------------------
#-------------------------------------------------------------------------------------------
LIMB = "without"
session = 'AUF'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_wSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_auf = pickle.load(file)
DATA_avg_auf = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_auf = pickle.load(file)
DATA_avg_pr_auf = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_auf = pickle.load(file)
session = 'ZU'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_wSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_wSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_wSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_wSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_zu = pickle.load(file)
DATA_avg_zu = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_zu = pickle.load(file)
DATA_avg_pr_zu = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_zu = pickle.load(file)
session = 'm1'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_woSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_woSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_woSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_woSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_m1 = pickle.load(file)
DATA_avg_m1 = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_m1 = pickle.load(file)
DATA_avg_pr_m1 = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_m1 = pickle.load(file)
session = 'm2'
file_path = f'../data/processed/matrices_data_{session}_{LIMB}_limbic.pkl'
data_avg_path = f'../data/processed/DATA_avg_woSC_{session}_{LIMB}_limb.csv'
data_sub_path = f'../data/processed/data_single_sub_woSC_{session}_{LIMB}_limb.pkl'
data_avg_pr_path = f'../data/processed/DATA_avg_woSC_pr_{session}_{LIMB}_limb.csv'
data_sub_pr_path = f'../data/processed/data_single_sub_woSC_pr_{session}_{LIMB}_limb.pkl'
with open(file_path, 'rb') as file:
data_input_m2 = pickle.load(file)
DATA_avg_m2 = pd.read_csv(data_avg_path)
with open(data_sub_path, 'rb') as file:
data_single_sub_m2 = pickle.load(file)
DATA_avg_pr_m2 = pd.read_csv(data_avg_pr_path)
with open(data_sub_pr_path, 'rb') as file:
data_single_sub_pr_m2 = pickle.load(file)
#--------------------------------------------------------------------------------------------
#------------------ concatenated all sessions data -------------------------------
#-------------------------------------------------------------------------------------------
df1 = zscore(data_single_sub_auf['ICallsub_w'], nan_policy='omit')
df2 = zscore(data_single_sub_zu['ICallsub_w'], nan_policy='omit')
df3 = zscore(data_single_sub_m1['ICallsub_w'], nan_policy='omit')
df4 = zscore(data_single_sub_m2['ICallsub_w'], nan_policy='omit')
ICallsub_allsessions = np.hstack((df1, df2, df3, df4))
#ICallsub_allsessions = df3
ICallsub_allsessions = pd.DataFrame(ICallsub_allsessions)
df1 = zscore(data_single_sub_pr_auf['degallsub_w'], nan_policy='omit')
df2 = zscore(data_single_sub_pr_zu['degallsub_w'], nan_policy='omit')
df3 = zscore(data_single_sub_pr_m1['degallsub_w'], nan_policy='omit')
df4 = zscore(data_single_sub_pr_m2['degallsub_w'], nan_policy='omit')
degallsub_allsessions = np.hstack((df1, df2, df3, df4))
#degallsub_allsessions = df3
degallsub_allsessions = pd.DataFrame(degallsub_allsessions)
df1 = zscore(data_input_auf['edata_medianallsub'], nan_policy='omit')
df2 = zscore(data_input_zu['edata_medianallsub'], nan_policy='omit')
df3 = zscore(data_input_m1['edata_medianallsub'], nan_policy='omit')
df4 = zscore(data_input_m2['edata_medianallsub'], nan_policy='omit')
edata_allsub_allsessions = np.hstack((df1, df2, df3, df4))
#edata_allsub_allsessions = df3
edata_allsub_allsessions = pd.DataFrame(edata_allsub_allsessions)
ICallsub_allsessions_avg = np.nanmean(ICallsub_allsessions, axis = 1)
degallsub_allsessions_avg = np.nanmean(degallsub_allsessions, axis = 1)
edata_allsub_allsessions_avg = np.nanmean(edata_allsub_allsessions, axis = 1)
#--------------------------------------------------------------------------------------------
#------------------ concatenated all sessions data -------------------------------
#-------------------------------------------------------------------------------------------
IC = ICallsub_allsessions_avg
DC = degallsub_allsessions_avg
CMRglc = edata_allsub_allsessions_avg
IC = np.array(IC).flatten()
DC = np.array(DC).flatten()
CMRglc = np.array(CMRglc).flatten()
alldata_avg = pd.DataFrame({'IC':IC, 'DC':DC, 'pet':CMRglc})
r_IC_CMRglc,p_IC_CMRglc = pcor(IC,CMRglc)
print('correlation MwC and CMRglc allsession datasets',r_IC_CMRglc )
r_DC_CMRglc , p_DC_CMRglc= pcor(DC,CMRglc)
print('correlation DC and CMRglc allsession datasets',r_DC_CMRglc)
r_IC_DC , p_IC_DC = pcor(IC,DC)
print('correlation MwC and DC allsession datasets',r_IC_DC)
sub_size = np.size(edata_allsub_allsessions, axis =1)
corr_IC_E_allsub = np.zeros((sub_size))
p_IC_E_allsub = np.zeros((sub_size))
for j in range(sub_size):
corr_IC_E_allsub[j], p_IC_E_allsub[j] = pcor(edata_allsub_allsessions.iloc[:,j], ICallsub_allsessions.iloc[:,j])
print('min p-val in subject level data:',p_IC_E_allsub.min())
print('max p-val in subject level data:',p_IC_E_allsub.max())
# # Load data
DC = degallsub_allsessions_avg
MwC = ICallsub_allsessions_avg
CMRglc = edata_allsub_allsessions_avg
# Flatten arrays
MwC = np.array(MwC).flatten()
DC = np.array(DC).flatten()
CMRglc = np.array(CMRglc).flatten()
# Create dataframe
df = pd.DataFrame({"MwC": MwC, "DC": DC, "CMRglc": CMRglc})
# Model 1: CMRglc ~ MCC (Baseline model)
X1 = sm.add_constant(df["MwC"]) # Add intercept
model1 = sm.OLS(df["CMRglc"], X1).fit()
# Model 2: CMRglc ~ MCC + DC (Full model)
X2 = sm.add_constant(df[["MwC", "DC"]]) # Add intercept
model2 = sm.OLS(df["CMRglc"], X2).fit()
# Print Model Summaries
print("\n### Model 1: CMRglc ~ MwC ###")
print(model1.summary())
print("\n### Model 2: CMRglc ~ MwC + DC ###")
print(model2.summary())
# Compare R² values
r2_model1 = model1.rsquared
r2_model2 = model2.rsquared
print(f"\nR² for Model 1 (MwC only): {r2_model1:.4f}")
print(f"R² for Model 2 (MwC + DC): {r2_model2:.4f}")
# ANOVA Comparison (F-test)
anova_results = sm.stats.anova_lm(model1, model2)
anova_p_value = anova_results["Pr(>F)"][1]
print(f"\nANOVA F-test p-value: {anova_p_value:.4f}")
# Interpretation
print(anova_results)
if anova_p_value < 0.01:
print("✅ Adding DC significantly improves the model. Both MwC and DC together explain more variance in CMRglc.")
else:
print("❌ Adding DC does not significantly improve the model over MwC alone.")
# %% [markdown]
# ## Steiger test DC vs MCC
# %%
n = 45
z_score, p_value = steiger_z_test(r_IC_CMRglc, r_DC_CMRglc, r_IC_DC, n)
print("steiger:", z_score, p_value)

Figures_codes.ipynb at commit 36f1937, under MIT · at the source

OverviewAuthors: Mahnaz Ashrafi1, Laura Fraticelli1,2, Gabriel Castrillón1,2,3, Valentin Riedl1,2
ORCID iDs: Mahnaz Ashrafi, Gabriel Castrillón, Valentin Riedl
Department of Neuroradiology, University Hospital TUM Klinikum, Technical University of Munich, Munich 81675, Germany
Department of Neuroradiology, Universitätsklinikum, Friedrich-Alexander-University Erlangen-Nuernberg, Erlangen 91054, Germany
Research Group in Medical Imaging, SURA Ayudas Diagnósticas, Medellín 6CQ4+M4, Colombia
Institutions: Technical University of Munich (Germany); Friedrich-Alexander-Universität Erlangen-Nürnberg (Germany)
Journal: Proceedings of the National Academy of Sciences of the United States of America, volume 123, issue 26, article e2531706123
Dates: received 4 November 2025; accepted 11 May 2026; published online 22 June 2026; in print 30 June 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1073/pnas.2531706123 · PMID 42330267 · PMCID PMC13321360 · OpenAlex W4414757148
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.10.02.679973
Status: code verified
Categories: PET / SPECT (modality), human (organism), other condition (population), cellular / molecular (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning, Preprocessing, Graphs, fMRI & imaging, Smoothing, state filtering, decompositions
Keywords: brain energy metabolism, brain connectome, synaptic integration, neurodegeneration
MeSH: Brain*, Connectome*, Neurodegenerative Diseases*, Synapses*, Energy Metabolism, Humans, Magnetic Resonance Imaging, Nerve Net, Positron-Emission Tomography (* major topic)
Journal subjects: Biological Sciences, Neuroscience
Topic: Mitochondrial Function and Pathology (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: European Research Council (759659)
Citations: cited by 1 paper (Europe PMC); 77 references in the paper
AbstractThe abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

SarahMorgan/Morphometric_Similarity_SZ

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 9134abfbc773c57a3a91b96612a7c8fa5969732b, 8 May 2019

Size: 7 files, 0 scripts

Software Heritage: not archived

Found in: “Data, Materials, and Software Availability”

Holds: README, license file

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

2 files

LICENSE — License, 21 lines
README.md — Text, 10 lines

At the source: github.com/SarahMorgan/Morphometric_Similarity_SZ

NeuroenergeticsLab/metabolism_weighted_centrality

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 36f1937bd610d762f1b92358c36158550f08e1d3, 8 June 2026

Languages: Python (5), Jupyter (3)

Size: 159 files, 8 scripts

Software Heritage: not archived

Found in: “Data, Materials, and Software Availability”

Holds: README, license file, environment (environment.yml, requirements.txt, setup.py), documentation, 3 notebooks

Not found: CITATION.cff, tests, continuous integration

Tools: Plotly (3 files), SciPy (3 files), BrainSMASH (2 files), igraph (2 files), Matplotlib (2 files), NetworkX (2 files), NiBabel (2 files), NumPy (2 files), pandas (2 files), statsmodels (2 files), Nilearn (1 file), scikit-learn (1 file), seaborn (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

10 files

docs/conf.py — Python, 244 lines
scripts/Figures_codes.ipynb — Jupyter, 4,803 lines, 5 matches
scripts/MwC_generation.ipynb — Jupyter, 166 lines, 1 match
setup.py — Python, 10 lines
src/Functions.ipynb — Jupyter, 1,727 lines, 3 matches
src/__init__.py — Python, 3 lines
src/functions.py — Python, 1,133 lines, 3 matches
test_environment.py — Python, 25 lines
LICENSE — License, 10 lines
README.md — Text, 50 lines

At the source: github.com/NeuroenergeticsLab/metabolism_weighted_centrality

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

Tracing map

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

What the map holds:

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

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

Data

Datasets cited

openneuro:ds004513 — at OpenNeuro; found in “Data, Materials, and Software Availability”

Code and data availability statement

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

it points to a dataset: OpenNeuro ds004513
it points to the authors' code: NeuroenergeticsLab/metabolism_weighted_centrality, SarahMorgan/Morphometric_Similarity_SZ

Read it in the paper: doi.org/10.1073/pnas.2531706123.

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, 4 authors, 4 keywords, 9 MeSH terms, 1 funder, 72 references.

Cite

This paper

Ashrafi, M., Fraticelli, L., Castrillón, G., & Riedl, V. (2026). Metabolism-weighted brain connectome reveals synaptic integration and vulnerability to neurodegeneration. Proceedings of the National Academy of Sciences of the United States of America, 123(26), e2531706123. https://doi.org/10.1073/pnas.2531706123

BibTeX

@article{ashrafi2026metabolism, author = {Ashrafi, Mahnaz and Fraticelli, Laura and Castrillón, Gabriel and Riedl, Valentin}, title = {{Metabolism-weighted brain connectome reveals synaptic integration and vulnerability to neurodegeneration}}, journal = {Proceedings of the National Academy of Sciences of the United States of America}, year = {2026}, month = jun, volume = {123}, number = {26}, pages = {e2531706123}, publisher = {National Academy of Sciences}, issn = {0027-8424}, doi = {10.1073/pnas.2531706123}, url = {https://doi.org/10.1073/pnas.2531706123}, pmid = {42330267}, pmcid = {PMC13321360} }

RIS

TY - JOUR AU - Ashrafi, Mahnaz AU - Fraticelli, Laura AU - Castrillón, Gabriel AU - Riedl, Valentin TI - Metabolism-weighted brain connectome reveals synaptic integration and vulnerability to neurodegeneration T2 - Proceedings of the National Academy of Sciences of the United States of America J2 - Proc Natl Acad Sci U S A PY - 2026 DA - 2026/06/22 VL - 123 IS - 26 SP - e2531706123 SN - 0027-8424 PB - National Academy of Sciences DO - 10.1073/pnas.2531706123 UR - https://doi.org/10.1073/pnas.2531706123 LA - en ER -

CSL-JSON

{ "id": "10.1073/pnas.2531706123", "type": "article-journal", "title": "Metabolism-weighted brain connectome reveals synaptic integration and vulnerability to neurodegeneration", "container-title": "Proceedings of the National Academy of Sciences of the United States of America", "author": [ { "family": "Ashrafi", "given": "Mahnaz" }, { "family": "Fraticelli", "given": "Laura" }, { "family": "Castrillón", "given": "Gabriel" }, { "family": "Riedl", "given": "Valentin" } ], "container-title-short": "Proc Natl Acad Sci U S A", "volume": "123", "issue": "26", "page": "e2531706123", "DOI": "10.1073/pnas.2531706123", "PMID": "42330267", "PMCID": "PMC13321360", "ISSN": "0027-8424", "publisher": "National Academy of Sciences", "URL": "https://doi.org/10.1073/pnas.2531706123", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 22 ] ] } }

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

Similar papers

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

[1] doi:10.1038/s42003-025-09444-3 [code]: Decoupling of neurophysiological activity from structure mirrors global microarchitectural and neuromodulatory trends.
Journal: Communications biology
In common: Nilearn, NiBabel, scikit-learn, 4 other tools, cellular / molecular, 11 references
[2] doi:10.1186/s12916-026-04903-y [code]: Structural connectome architecture and biological vulnerability shape cortical atrophy in cocaine use disorder.
Journal: BMC medicine
In common: Nilearn, statsmodels, seaborn, 5 other tools, other condition, 9 references
[3] doi:10.1038/s41398-026-04025-2 [code]: Brain energetic landscapes shape state dysregulation in major depressive disorder: a morphological network controllability perspective.
Journal: Translational psychiatry
In common: BrainSMASH, Nilearn, NiBabel, 6 other tools, 6 references
[4] doi:10.1162/imag.a.1275 [code]: Glucose metabolism echoes long-range temporal correlations in the human brain.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: BrainSMASH, NiBabel, scikit-learn, 4 other tools, PET / SPECT, 6 references
[5] doi:10.1038/s41467-026-75959-w [code]: Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: NetworkX, Nilearn, statsmodels, 7 other tools, 6 references
[6] doi:10.1038/s41467-026-71270-w [code]: Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan.
Journal: Nature communications
In common: NetworkX, Nilearn, NiBabel, 6 other tools, 6 references
[7] doi:10.1186/s12916-026-04978-7 [code]: Mapping shared and specific cortical after-effects of repetitive TMS on brain function.
Journal: BMC medicine
In common: BrainSMASH, Nilearn, NiBabel, 5 other tools, 5 references
[8] doi:10.1038/s41398-026-03965-z [code]: Disentangling individual heterogeneity reveals robust network and molecular signatures of major depressive disorder with suicidal ideation.
Journal: Translational psychiatry
In common: NetworkX, statsmodels, NiBabel, 6 other tools, cellular / molecular, 5 references
[9] doi:10.7554/elife.103097 [code]: Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: Nilearn, Plotly, statsmodels, 5 other tools, 6 references
[10] doi:10.1038/s41467-026-76812-w [code]: Assessing molecular, cellular and transcriptomic bases of laminar perfusion and cytoarchitecture coupling in the human cortex.
Journal: Nature communications
In common: Nilearn, NiBabel, scikit-learn, 4 other tools, cellular / molecular, 6 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.