A functional connectome identifies vulnerability to cognitive decline and Alzheimer’s disease progression

Code ↔ Paper

3 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 3 matches

[1] § MATERIALS AND METHODS › Neuropsychological assessments ↔ clustering.py, lines 1478 to 1533 · score 0.90 · CFT VEG, RAVLT Immediate, RAVLT Learning, Forgetting, bills, events · Report a mismatch
[2] § MATERIALS AND METHODS › Classification between SA and AD ↔ SA_AD_classification.py, lines 568 to 628 · score 0.58 · AD classification, L1, penalty, probabilities, AUC, accuracy · Report a mismatch
[3] § MATERIALS AND METHODS › Superager definition ↔ clustering.py, lines 1478 to 1533 · score 0.55 · CFT VEG, vegetable, delayed, TMT, RAVLT, animal · Report a mismatch

Paper

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

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,691 lines · 97 KB · no license · 2 matches

# -*- coding: utf-8 -*-
"""
Created on Mon Feb 19 00:07:12 2024
@author: 99488
"""
'''
result with age corrected subjects from adni, then cross validation in subject level were obtained, combat all subjects
'''
import os
from arc_harmonization import FoldCombat, load_adni_calibration, participant_scores
import pandas as pd
import numpy as np
import scipy.io as sio
from pathlib import Path
from bids.layout import parse_file_entities
import sys
sys.path.append(r'F:\OPT\research\fMRI\utils_for_all')
from common_utils import read_singal_fmri_ts, get_fc, t_test, keep_triangle_half, sens_spec, setup_seed, \
vector_to_matrix, heatmap, get_net_net_connect, diverge_map_unipolar_colore, cohen_d
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import StratifiedKFold
from sklearn.utils.class_weight import compute_sample_weight
from sklearn.svm import LinearSVC, SVC, OneClassSVM
from sklearn.metrics import accuracy_score, r2_score, auc, roc_curve, roc_auc_score
from sklearn.ensemble import RandomForestClassifier, IsolationForest
import xgboost as xgb
from mrmr import mrmr_classif
import matplotlib.pyplot as plt
from sklearn.impute import SimpleImputer
import copy
from scipy.stats import pearsonr, chi2_contingency, shapiro, kstest, iqr, ranksums, kendalltau, f_oneway, kruskal, kendalltau, spearmanr
import matplotlib
import seaborn as sns
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
from collections import namedtuple, OrderedDict
from matplotlib import colors as mcolors
from statsmodels.stats.multitest import multipletests
import statsmodels.api as sm
from statsmodels.stats.anova import anova_lm
from sklearn.linear_model import LinearRegression
from statsmodels.formula.api import ols
font2 = {'family': 'Tahoma', 'weight': 'bold', 'size': 30}
matplotlib.rc('font', **font2)
setup_seed(6)
Color = namedtuple('RGB', 'red, green, blue')
colors = {} # dict of colors
class RGB(Color):
def hex_format(self):
return '#{:02X}{:02X}{:02X}'.format(self.red, self.green, self.blue)
# Color Contants
red = RGB(0.941, 0.502, 0.498)
blue = RGB(0.388, 0.584, 0.933)
def make_colormap(seq):
"""Return a LinearSegmentedColormap
seq: a sequence of floats and RGB-tuples. The floats should be increasing
and in the interval (0,1).
"""
seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3]
cdict = {'red': [], 'green': [], 'blue': []}
for i, item in enumerate(seq):
if isinstance(item, float):
r1, g1, b1 = seq[i - 1]
r2, g2, b2 = seq[i + 1]
cdict['red'].append([item, r1, r2])
cdict['green'].append([item, g1, g2])
cdict['blue'].append([item, b1, b2])
return mcolors.LinearSegmentedColormap('CustomMap', cdict)
def diverge_map(high=(0.565, 0.392, 0.173), low=(0.094, 0.310, 0.635)):
'''
low and high are colors that will be used for the two
ends of the spectrum. they can be either color strings
or rgb color tuples
'''
c = mcolors.ColorConverter().to_rgb
if isinstance(low, str): low = c(low)
if isinstance(high, str): high = c(high)
return make_colormap([low, c('white'), 0.5, c('white'), high])
# ---- Load cohort inputs (AD=1, SA=0, CN=2, MCI=3) ----
info = sio.loadmat(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\superage_data_adni.mat')
sa_dem_fc_adni = info['fc_sa_dem']
sa_dem_dx_adni = np.squeeze(info['dx_sa_dem'])
sa_dem_age_adni = np.squeeze(info['age_sa_dem'])
sa_dem_sex_adni = np.squeeze(info['sex_sa_dem'])
sa_dem_sex_adni = np.array([2 if sex == 'Female' else 1 for sex in sa_dem_sex_adni])
sa_dem_sub_adni = np.squeeze(info['sub_id_sa_dem'])
sa_dem_sess_adni = np.squeeze(info['sess_sa_dem'])
df=pd.DataFrame(np.c_[sa_dem_sub_adni, sa_dem_sess_adni],columns=['sub','sess'])
df['sub_sess']=df['sub']+'_'+df['sess']
sa_dem_sub_sess_adni = df['sub_sess'].values
half_fc_adni = keep_triangle_half(sa_dem_fc_adni.shape[1] * (sa_dem_fc_adni.shape[1]-1)//2, sa_dem_fc_adni.shape[0], sa_dem_fc_adni)
all_fc_adni = info['adni_fc_all']
all_dx_adni = np.squeeze(info['adni_dx_all'])
all_age_adni = np.squeeze(info['adni_age_all'])
all_sex_adni = np.squeeze(info['adni_sex_all'])
# all_sex_adni = np.array([2 if sex == 'Female' else 1 for sex in all_sex_adni])
all_sub_adni = np.squeeze(info['adni_sub_all'])
all_sess_adni = np.squeeze(info['adni_sess_all'])
other_adni_sub = np.setdiff1d(all_sub_adni, sa_dem_sub_adni)
other_fc_adni = []
other_dx_adni = []
other_age_adni = []
other_sex_adni = []
other_sess_adni = []
other_sub_adni = []
for i,sub in enumerate(other_adni_sub):
other_fc_adni.extend(list(all_fc_adni[all_sub_adni==sub]))
other_dx_adni.extend(list(all_dx_adni[all_sub_adni==sub]))
other_age_adni.extend(list(all_age_adni[all_sub_adni==sub]))
other_sex_adni.extend(list(all_sex_adni[all_sub_adni==sub]))
other_sess_adni.extend(list(all_sess_adni[all_sub_adni==sub]))
other_sub_adni.extend(list(all_sub_adni[all_sub_adni==sub]))
other_fc_adni = np.array(other_fc_adni)
other_dx_adni = np.array(other_dx_adni)
other_age_adni = np.array(other_age_adni)
other_sex_adni = np.array(other_sex_adni)
other_sess_adni = np.array(other_sess_adni)
other_sub_adni = np.array(other_sub_adni)
low_qua_fc_idx = np.where(np.isnan(np.sum(np.sum(other_fc_adni,-1),-1)))[0]
other_fc_adni = np.delete(other_fc_adni, low_qua_fc_idx, 0)
other_dx_adni = np.delete(other_dx_adni, low_qua_fc_idx, 0)
other_age_adni = np.delete(other_age_adni, low_qua_fc_idx, 0)
other_sex_adni = np.delete(other_sex_adni, low_qua_fc_idx, 0)
other_sess_adni = np.delete(other_sess_adni, low_qua_fc_idx, 0)
other_sub_adni = np.delete(other_sub_adni, low_qua_fc_idx, 0)
df=pd.DataFrame(np.c_[other_sub_adni, other_sess_adni],columns=['sub','sess'])
df['sub_sess']=df['sub']+'_'+df['sess']
other_sub_sess_adni = df['sub_sess'].values
for i in range(len(other_sex_adni)):
if 'F' in other_sex_adni[i]:
other_sex_adni[i] = 2
elif 'M' in other_sex_adni[i]:
other_sex_adni[i] = 1
if 'CN' in other_dx_adni[i]:
other_dx_adni[i] = 2
elif 'MCI' in other_dx_adni[i]:
other_dx_adni[i] = 3
elif 'Dementia' in other_dx_adni[i]:
other_dx_adni[i] = 1
other_dx_adni = other_dx_adni.astype(float)
other_sex_adni = other_sex_adni.astype(float)
info_oas = sio.loadmat(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\superage_data_oas_dem.mat')
sa_dem_fc_oas = info_oas['multi_run_fcs_sa_dem']
sa_dem_dx_oas = np.squeeze(info_oas['multi_run_dx_sa_dem'])
sa_dem_sub_oas_ = np.squeeze(info_oas['multi_run_subs_sa_dem'])
sa_dem_sess_oas = np.squeeze(info_oas['multi_run_sess_sa_dem'])
all_pcd_oas = np.squeeze(info_oas['multi_run_pcd_all'])
sa_dem_sub_oas = pd.Series(sa_dem_sub_oas_).str.cat(pd.Series(sa_dem_sess_oas), sep = '_').values
sex_oas = np.squeeze(info_oas['sex_sa_dem'])
pcd_oas = info_oas['pcd_sa_dem']
sess_oas = info_oas['sess_sa_dem']
age_oas = pcd_oas[:,2].astype(float).squeeze()
sub_oas = pcd_oas[:,0]
sub_oas = np.array([sub[0].split('_')[0] for sub in sub_oas])
oasis_path = r'E:\PHD\learning\research\AD\data\OASIS3'
oasis_demo = pd.read_csv(os.path.join(oasis_path, 'UDS-A2-InformatDemos.csv'))
oasis_demo2 = pd.read_csv(os.path.join(oasis_path, 'OASIS3_UDSa2_cs_demo.csv'))
oasis_demo_sub = oasis_demo['Subject']
oasis_demo_sub2 = oasis_demo2['OASISID']
sa_dem_sex_oas = []
sa_dem_age_oas = []
for i, sub in enumerate(sa_dem_sub_oas_):
# sa_dem_sex_oas.append(sex_oas[(sub_oas==sub)&(sess_oas == sa_dem_sess_oas[i])])
if sum(oasis_demo_sub == sub) >= 1:
sa_dem_sex_oas.append(oasis_demo['INSEX'][oasis_demo_sub == sub].values[0])
else:
if sum(oasis_demo_sub2 == sub) >= 1:
sa_dem_sex_oas.append(oasis_demo2['INSEX'][oasis_demo_sub2 == sub].values[0])
sa_dem_age_oas.append(age_oas[(sub_oas==sub)&(sess_oas == sa_dem_sess_oas[i])])
sa_dem_sex_oas = np.array(sa_dem_sex_oas).squeeze() # female: 2; male: 1
sa_dem_age_oas = np.array(sa_dem_age_oas).squeeze()
all_fc_oas = info_oas['multi_run_fcs_all']
all_sub_oas = np.squeeze(info_oas['multi_run_subs_all'])
all_sess_oas = np.squeeze(info_oas['multi_run_sess_all'])
all_pcd_oas = np.squeeze(info_oas['multi_run_pcd_all'])
all_sex_oas = np.squeeze(info_oas['multi_run_sex_all'])
all_age_oas = all_pcd_oas[:,2].astype(float).squeeze()
all_dx_oas = all_pcd_oas[:,1]
other_oas_sub = np.setdiff1d(all_sub_oas, sa_dem_sub_oas_)
other_fc_oas = []
other_dx_oas = []
other_age_oas = []
other_sex_oas = []
other_sess_oas = []
other_sub_oas_multises = []
for i,sub in enumerate(other_oas_sub):
other_fc_oas.extend(list(all_fc_oas[all_sub_oas==sub]))
other_dx_oas.extend([s[0][0] if isinstance(s[0], np.ndarray) else s[0] for s in all_dx_oas[all_sub_oas==sub]])
other_age_oas.extend(list(all_age_oas[all_sub_oas==sub]))
other_sex_oas.extend(list(all_sex_oas[all_sub_oas==sub]))
other_sess_oas.extend(list(all_sess_oas[all_sub_oas==sub]))
other_sub_oas_multises.extend(list(all_sub_oas[all_sub_oas==sub]))
other_fc_oas = np.array(other_fc_oas)
other_dx_oas = np.array(other_dx_oas)
other_age_oas = np.array(other_age_oas)
other_sex_oas = np.array(other_sex_oas)
other_sess_oas = np.array(other_sess_oas)
other_sub_oas_multises = np.array(other_sub_oas_multises)
mask_nc = pd.Series(other_dx_oas).str.contains('normal', regex=False)
mask_ad = (pd.Series(other_dx_oas).str.contains('AD', regex=False)) |(pd.Series(other_dx_oas).str.contains('DAT', regex=False)) | \
(pd.Series(other_dx_oas)=='Frontotemporal demt. prim') | (pd.Series(other_dx_oas)=='DLBD, primary') | (pd.Series(other_dx_oas).str.contains('primary', regex=False)) | \
(pd.Series(other_dx_oas)=='Vascular Demt, secondary') | (pd.Series(other_dx_oas)=='uncertain dementia') | \
(pd.Series(other_dx_oas)=='uncertain, possible NON AD dem') | (pd.Series(other_dx_oas)=='DLBD, primary')
mask_mci = pd.Series(other_dx_oas).str.contains('0.5 in memory', regex=False) | (pd.Series(other_dx_oas)=='Unc: ques. Impairment')
other_dx_oas_encode = np.zeros((len(other_dx_oas)))
other_dx_oas_encode[mask_nc] = 2
other_dx_oas_encode[mask_mci] = 3
other_dx_oas_encode[mask_ad] = 1
low_qua_fc_idx = np.where(np.isnan(np.sum(np.sum(other_fc_oas,-1),-1)))[0]
other_fc_oas = np.delete(other_fc_oas, low_qua_fc_idx, 0)
other_dx_oas_encode_raw = np.delete(other_dx_oas_encode, low_qua_fc_idx, 0)
other_age_oas_raw = np.delete(other_age_oas, low_qua_fc_idx, 0)
other_sex_oas_raw = np.delete(other_sex_oas, low_qua_fc_idx, 0)
other_sess_oas_raw = np.delete(other_sess_oas, low_qua_fc_idx, 0)
other_sub_oas_multises = np.delete(other_sub_oas_multises, low_qua_fc_idx, 0)
half_fc_oas_sa_dem = keep_triangle_half(sa_dem_fc_oas.shape[1] * (sa_dem_fc_oas.shape[1]-1)//2, sa_dem_fc_oas.shape[0], sa_dem_fc_oas)
half_fc_oas_other = keep_triangle_half(other_fc_oas.shape[1] * (other_fc_oas.shape[1]-1)//2, other_fc_oas.shape[0], other_fc_oas)
half_fc_adni_other = keep_triangle_half(other_fc_adni.shape[1] * (other_fc_adni.shape[1]-1)//2, other_fc_adni.shape[0], other_fc_adni)
info_habs = sio.loadmat(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\superage_data_habs.mat')
sa_pcd_habs = info_habs['multi_run_pcd_sa'].squeeze()
sa_fc_habs = info_habs['multi_run_fcs_sa'][:,:100,:100]
sa_sub_sess_habs = np.array([sub[0] for sub in sa_pcd_habs[:,-1].squeeze()])
sa_sex_habs = np.array([2 if sex[0] == 'F' else 1 for sex in sa_pcd_habs[:,59]])
sa_age_habs = np.array([sub[0][0] for sub in sa_pcd_habs[:,58].squeeze()])
sa_sex_habs = sa_sex_habs.astype(float)
sa_dx_habs = np.zeros((len(sa_sex_habs)))
half_fc_habs = keep_triangle_half(sa_fc_habs.shape[1] * (sa_fc_habs.shape[1]-1)//2, sa_fc_habs.shape[0], sa_fc_habs)
sa_sub_habs = pd.Series(sa_sub_sess_habs).str.split('_', expand = True)[0]
habs_pcd_all = info_habs['multi_run_pcd_all'].squeeze()
habs_fc_all = info_habs['multi_run_fc_all'][:,:100,:100]
habs_sub_sess_all = np.array([sub[0] for sub in habs_pcd_all[:,-1].squeeze()])
habs_sex_all = np.array([2 if sex[0] == 'F' else 1 for sex in habs_pcd_all[:,59]])
habs_age_all = np.array([sub[0][0] for sub in habs_pcd_all[:,58].squeeze()])
habs_sex_all = habs_sex_all.astype(float)
habs_dx_all = np.array([sub[0] for sub in habs_pcd_all[:,67].squeeze()])
habs_dx_all[habs_dx_all == 'CN'] = 2
habs_dx_all[habs_dx_all == 'MCI'] = 3
habs_dx_all[habs_dx_all == 'Dementia'] = 1
habs_dx_all = habs_dx_all.astype(float)
half_habs_fc_all = keep_triangle_half(habs_fc_all.shape[1] * (habs_fc_all.shape[1]-1)//2, habs_fc_all.shape[0], habs_fc_all)
habs_sub_all = pd.Series(habs_sub_sess_all).str.split('_', expand = True)[0]
other_habs_sub = np.unique(habs_sub_sess_all[~np.isin(np.asarray(habs_sub_all), np.asarray(sa_sub_habs))])
other_fc_habs = []
other_dx_habs = []
other_age_habs = []
other_sex_habs = []
other_sub_sess_habs = []
for i,sub in enumerate(other_habs_sub):
other_fc_habs.extend(list(half_habs_fc_all[habs_sub_sess_all==sub]))
other_dx_habs.extend(list(habs_dx_all[habs_sub_sess_all==sub]))
other_age_habs.extend(list(habs_age_all[habs_sub_sess_all==sub]))
other_sex_habs.extend(list(habs_sex_all[habs_sub_sess_all==sub]))
other_sub_sess_habs.extend(list(habs_sub_sess_all[habs_sub_sess_all==sub]))
other_fc_habs = np.array(other_fc_habs)
other_dx_habs = np.array(other_dx_habs)
other_age_habs = np.array(other_age_habs)
other_sex_habs = np.array(other_sex_habs)
other_sub_sess_habs = np.array(other_sub_sess_habs)
half_fc_ind = np.r_[half_fc_oas_sa_dem, half_fc_habs]
age_ind = np.r_[sa_dem_age_oas, sa_age_habs]
sex_ind = np.r_[sa_dem_sex_oas, sa_sex_habs]
dx_ind = np.r_[sa_dem_dx_oas, sa_dx_habs]
sub_ind = np.r_[sa_dem_sub_oas, sa_sub_sess_habs]
sub_ind_ = np.r_[sa_dem_sub_oas_, sa_sub_habs]
site_ind = np.r_[np.zeros(len(sa_dem_sub_oas)), sa_dx_habs+20]
# sub_ind is visit-level; the original subject IDs live in sub_ind_.
participant_keys = np.array([f'{site}:{str(pid).strip()}' for site, pid in zip(site_ind, sub_ind_)])
subject_used, uni_id = np.unique(participant_keys, return_index=True)
if any(np.unique(dx_ind[participant_keys == person]).size != 1 for person in subject_used):
raise ValueError('Anchor diagnosis changes across visits; verify anchor inclusion first')
print('Counts of raw subjects: {}, superager: {}, dementia: {}'.format(len(set(sub_ind)), sum(dx_ind[uni_id]==0), sum(dx_ind[uni_id]==1)))
print('Counts of used runs: {}, superager: {}, dementia: {}'.format(len(sub_ind), sum(dx_ind==0), sum(dx_ind==1)))
scaler = StandardScaler()
half_fc_adni = scaler.fit_transform(half_fc_adni.T).T
half_fc_ind = scaler.fit_transform(half_fc_ind.T).T
half_fc_oas_other = scaler.fit_transform(half_fc_oas_other.T).T
half_fc_adni_other = scaler.fit_transform(half_fc_adni_other.T).T
other_fc_habs = scaler.fit_transform(other_fc_habs.T).T
df=pd.DataFrame(np.c_[other_sub_oas_multises, other_sess_oas_raw],columns=['sub','sess'])
df['sub_sess']=df['sub']+'_'+df['sess']
sub_sess_oas = df['sub_sess'].values
sub_combine_uni_oas, uni_id_oas = np.unique(df['sub_sess'], return_index=True)
#################averaging runs
half_fc_oas_other_uni = []
for i, sub in enumerate(sub_combine_uni_oas):
mask = df['sub_sess'].values == sub
half_fc_oas_other_uni.append(half_fc_oas_other[mask].mean(0))
half_fc_oas_other_uni = np.array(half_fc_oas_other_uni)
other_dx_oas_encode = other_dx_oas_encode_raw[uni_id_oas]
other_age_oas = other_age_oas_raw[uni_id_oas]
other_sex_oas = other_sex_oas_raw[uni_id_oas]
#####################redefine MCI based on ADNI criteria
oasis_sa_sub_sess_uniq, oasis_sa_sub_sess_idx, oasis_sa_sub_sess_idx_reverse = np.unique(sub_combine_uni_oas, return_index=True, return_inverse=True)
oas_sess_info = sio.loadmat(r'F:\PHD\learning\project\super_age\sess_match_oas_all.mat')
day_all = oas_sess_info['day']
sess_all = oas_sess_info['sess']
idx = np.argsort(sess_all)
sess_all = sess_all[idx]
day_all = pd.Series(day_all[idx]).str.split('d', expand = True).iloc[:,-1].astype(int)
oas_pcd_info = sio.loadmat(r'E:\PHD\learning\research\AD\data\OASIS3\fMRIdata_OASIS3_struct.mat')
tvar = oas_pcd_info['OASIS3_Phenotypic_struct']
PCD_data = []
pcd_keys = ['subjectID_Date', 'apoe', 'NPIQINF', 'NPIQINFX', 'DEL',
'DELSEV', 'HALL', 'HALLSEV', 'AGIT', 'AGITSEV', 'DEPD', 'DEPDSEV', 'ANX',
'ANXSEV', 'ELAT', 'ELATSEV', 'APA', 'APASEV', 'DISN', 'DISNSEV', 'IRR',
'IRRSEV', 'MOT', 'MOTSEV', 'NITE', 'NITESEV', 'APP', 'APPSEV', 'mmse',
'cdr', 'commun','homehobb', 'judgment', 'memory', 'orient', 'perscare', 'sumbox',
'DIGIF', 'DIGIB', 'ANIMALS', 'VEG', 'TRAILA', 'TRAILALI', 'TRAILB',
'TRAILBLI', 'WAIS', 'LOGIMEM', 'MEMUNITS', 'MEMTIME', 'BOSTON']
for idx in range(len(tvar)):
sub_data = []
for name in pcd_keys:
value = tvar[idx][0][name]
for i in range(3):
if isinstance(value, np.ndarray):
try:
value = value[0]
except IndexError:
if name == 'DX':
value = 'nan'
else:
value = "error"
else:
if name == 'subjectID_Date' or name == 'dx1':
value = str(value)
sub_data.append(value)
break
PCD_data.append(sub_data)
PCD_data = np.array(PCD_data)
PCD_data[:,1] = pd.Series(PCD_data[:,1]).str.count('4')
#update NPI 0 score
for i in [5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27]:
mask = PCD_data[:, i] == 'nan'
PCD_data[np.isfinite(PCD_data[:,2].astype(float))&mask, i] = 0
PCD_data[np.isfinite(PCD_data[:,2].astype(float))][:,[5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27]] = 0
PCD_sub = np.array([sub[:8] for sub in PCD_data[:,0]])
PCD_day = np.array([sub[10:] for sub in PCD_data[:,0]]).astype(int)
PCD_sess = []
day_min = []
for da in PCD_day:
mask = day_all==da
if sum(mask) > 0:
PCD_sess.append(sess_all[day_all==da][0])
else:
diff = day_all-da
diff_min_id = np.argmin(abs(diff))
diff_min = diff[diff_min_id]
day_min.append(diff_min)
PCD_sess.append(sess_all[diff_min_id])
PCD_sess = np.array(PCD_sess)
PCD_interest = []
for sub_sess in oasis_sa_sub_sess_uniq:
sub = sub_sess[:8]
sess = sub_sess[9:]
if sum((PCD_sub == sub) & (PCD_sess == sess))== 1:
PCD_interest.append(PCD_data[:,np.r_[[1, 2, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25], np.arange(27,50)]][
(PCD_sub == sub) & (PCD_sess == sess)].squeeze())
elif sum((PCD_sub == sub) & (PCD_sess == sess))> 1:
pcd = PCD_data[:,np.r_[[1, 2, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25], np.arange(27,50)]][
(PCD_sub == sub) & (PCD_sess == sess)].squeeze()
pcd_new = []
for i in range(pcd.shape[-1]):
info = pcd[:,i]
if np.isnan(info.astype(float)).sum() == len(info):
pcd_new.append(np.nan)
else:
pcd_new.append(info[~np.isnan(info.astype(float))][0].squeeze())
PCD_interest.append(np.array(pcd_new))
else:
PCD_interest.append(np.array([np.nan for i in range(len(['apoe', 'NPIQINF', 'DELSEV', 'HALLSEV', 'AGITSEV', 'DEPDSEV',
'ANXSEV', 'ELATSEV', 'APASEV', 'DISNSEV', 'IRRSEV', 'MOTSEV', 'NITESEV',
'APPSEV', 'mmse', 'cdr', 'commun','homehobb', 'judgment', 'memory', 'orient', 'perscare', 'sumbox',
'DIGIF', 'DIGIB', 'ANIMALS', 'VEG', 'TRAILA', 'TRAILALI', 'TRAILB',
'TRAILBLI', 'WAIS', 'LOGIMEM', 'MEMUNITS', 'MEMTIME', 'BOSTON']))]))
PCD_interest = np.array(PCD_interest)
df=pd.DataFrame(np.c_[other_sub_adni, other_sess_adni],columns=['sub','sess'])
df['sub_sess']=df['sub']+'_'+df['sess']
sub_combine_uni, uni_id = np.unique(df['sub_sess'], return_index=True)
half_fc_adni_other_uni = []
for i, sub in enumerate(sub_combine_uni):
mask = df['sub_sess'].values == sub
half_fc_adni_other_uni.append(half_fc_adni_other[mask].mean(0))
half_fc_adni_other_uni = np.array(half_fc_adni_other_uni)
other_dx_adni = other_dx_adni[uni_id]
other_sub_sess_adni = other_sub_sess_adni[uni_id]
other_age_adni = other_age_adni[uni_id]
other_age_adni[other_age_adni == 0] = np.nan
other_sex_adni = other_sex_adni[uni_id]
other_sub_adni = other_sub_adni[uni_id]
other_sess_adni = other_sess_adni[uni_id]
# Inputs to fold-local ComBat, after the original per-scan normalization.
half_fc_adni_correct = half_fc_adni
half_fc_adni_other_correct = half_fc_adni_other_uni
half_fc_oas_other_correct = half_fc_oas_other_uni
half_fc_habs_other_correct = other_fc_habs
half_fc_adni_correct_all = np.r_[half_fc_adni_correct, half_fc_adni_other_correct]
half_fc_adni_dx_all = np.r_[sa_dem_dx_adni, other_dx_adni]
half_fc_adni_age_all = np.r_[sa_dem_age_adni, other_age_adni]
half_fc_adni_sex_all = np.r_[sa_dem_sex_adni, other_sex_adni]
half_fc_adni_sub_sess_all = np.r_[sa_dem_sub_sess_adni, other_sub_sess_adni]
info_other_sa_define = sio.loadmat(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\superage_data_adni_logicmem_based.mat')
sa_dem_sub_adni_new_define = np.squeeze(info_other_sa_define['sub_id_sa_dem'])
sa_dem_sess_adni_new_define = np.squeeze(info_other_sa_define['sess_sa_dem'])
sa_dem_dx_adni_new_define = np.squeeze(info_other_sa_define['dx_sa_dem'])
df=pd.DataFrame(np.c_[sa_dem_sub_adni_new_define, sa_dem_sess_adni_new_define],columns=['sub','sess'])
df['sub_sess']=df['sub']+'_'+df['sess']
sa_dem_sub_sess_adni_new_define = df['sub_sess'].values
_, sa_adni_new_define_idx, sa_adni_idx = np.intersect1d(sa_dem_sub_sess_adni_new_define[sa_dem_dx_adni_new_define == 0], half_fc_adni_sub_sess_all[half_fc_adni_dx_all == 0], return_indices=True)
sa_fc_adni_new_define = half_fc_adni_correct_all[half_fc_adni_dx_all == 0][sa_adni_idx]
sa_sex_adni_new_define = half_fc_adni_sex_all[half_fc_adni_dx_all == 0][sa_adni_idx]
sa_age_adni_new_define = half_fc_adni_age_all[half_fc_adni_dx_all == 0][sa_adni_idx]
sa_sub_sess_adni_new_define = half_fc_adni_sub_sess_all[half_fc_adni_dx_all == 0][sa_adni_idx]
dem_fc_adni_new_define = half_fc_adni_correct[sa_dem_dx_adni == 1]
dem_sex_adni_new_define = sa_dem_sex_adni[sa_dem_dx_adni == 1]
dem_age_adni_new_define = sa_dem_age_adni[sa_dem_dx_adni == 1]
dem_sub_sess_adni_new_define = sa_dem_sub_sess_adni[sa_dem_dx_adni == 1]
sa_dem_fc_adni_new_define = np.r_[dem_fc_adni_new_define, sa_fc_adni_new_define]
sa_dem_sex_adni_new_define = np.r_[dem_sex_adni_new_define, sa_sex_adni_new_define]
sa_dem_age_adni_new_define = np.r_[dem_age_adni_new_define, sa_age_adni_new_define]
sa_dem_sub_sess_adni_new_define = np.r_[dem_sub_sess_adni_new_define, sa_sub_sess_adni_new_define]
sa_dem_sub_adni_new_define = pd.Series(sa_dem_sub_sess_adni_new_define).str.split(' ', expand = True).iloc[:,0]
sa_dem_dx_adni_new_define = np.r_[np.ones(len(dem_fc_adni_new_define)), np.zeros(len(sa_fc_adni_new_define))]
########################
# sub_ind is visit-level; the original subject IDs live in sub_ind_.
participant_keys = np.array([f'{site}:{str(pid).strip()}' for site, pid in zip(site_ind, sub_ind_)])
subject_used, uni_id = np.unique(participant_keys, return_index=True)
if any(np.unique(dx_ind[participant_keys == person]).size != 1 for person in subject_used):
raise ValueError('Anchor diagnosis changes across visits; verify anchor inclusion first')
dx_ind_uniq = dx_ind[uni_id]
site_ind_uniq = site_ind[uni_id]
site_ind_uniq[site_ind_uniq == 1] = 20
# half_fc_adni_correct = copy.deepcopy(half_fc_adni)
t_p_vector = np.zeros((4950, 2))
for i in range(4950):
t_p_vector[i, 0], t_p_vector[i, 1] = t_test(half_fc_ind[dx_ind == 1, i], half_fc_ind[dx_ind == 0, i], 0.05, False)
weight_sys, _ = vector_to_matrix(t_p_vector[:,0])
label_idx = np.arange(0,100,1)
cbarlabel = 'Weights'
plt.figure(figsize =(15,15))
ax = plt.gca()
im, cbar = heatmap(weight_sys.squeeze(), label_idx, ax=ax, cmap='RdBu_r', connect_type='roi',
cbarlabel=cbarlabel, half_or_full = 'half', with_diag = True)
folds = 10
kf = StratifiedKFold(n_splits=folds, random_state=9, shuffle=True)
perform = []
perform_sub = []
perform_oas = []
perform_habs = []
weights = []
t_p_vector_all = []
indep_pred_prob = []
indep_performance_all = []
indep_performance_logic = []
indep_pred_prob_logic = []
other_adni_pred_label = []
other_oas_pred_label = []
other_habs_pred_label = []
# Auxiliary ADNI contrasts; descriptive, not a new fitted classifier.
cn = half_fc_adni_other_correct[other_dx_adni == 2]
ad = half_fc_adni_correct[sa_dem_dx_adni == 1]
sa = half_fc_adni_correct[sa_dem_dx_adni == 0]
half_fc_cn_dem_pair = np.r_[cn, ad]
dx_cn_dem_pair = np.r_[np.zeros(len(cn)), np.ones(len(ad))]
half_fc_cn_sa_pair = np.r_[cn, sa]
dx_cn_sa_pair = np.r_[np.ones(len(cn)), np.zeros(len(sa))]
indep_cn_dem_prob = []
indep_cn_sa_prob = []
# ---- Participant CV, fold harmonization, and external scoring ----
adni_calibration = load_adni_calibration(np.r_[all_sub_adni, sa_dem_sub_adni], half_fc_ind.shape[1])
for train_index, test_index in kf.split(dx_ind_uniq, dx_ind_uniq):
tr_all_id = []
for sub in subject_used[train_index]:
tr_all_id.extend(list(np.where(participant_keys==sub)[0]))
tr_all_id = np.array(tr_all_id)
te_all_id = []
for sub in subject_used[test_index]:
te_all_id.extend(list(np.where(participant_keys==sub)[0]))
te_all_id = np.array(te_all_id)
if set(participant_keys[tr_all_id]) & set(participant_keys[te_all_id]):
raise RuntimeError('Participant appears in train and test')
te_subs = participant_keys[te_all_id]
te_sub_unqiu, te_id = np.unique(te_subs, return_index=True)
X_train = half_fc_ind[tr_all_id]
X_test = half_fc_ind[te_all_id]
y_train = dx_ind[tr_all_id]
y_test = dx_ind[te_all_id]
site_test = site_ind[te_all_id]
site_train = site_ind[tr_all_id]
site_test_te = site_test[te_id]
idx = np.arange(len(X_train))
np.random.default_rng(9 + len(weights)).shuffle(idx)
site_train = site_train[idx]
train_people = participant_keys[tr_all_id][idx]
X_train = X_train[idx]
y_train = y_train[idx]
p_thre = 0.01
people, inverse = np.unique(train_people, return_inverse=True)
subject_x = np.stack([X_train[inverse == j].mean(axis=0) for j in range(len(people))])
subject_y = np.array([y_train[inverse == j][0] for j in range(len(people))])
person_sites = np.array([site_train[inverse == j][0] for j in range(len(people))])
# The separate ADNI calibration set defines batch 2 in the shared space.
combat = FoldCombat().fit(np.r_[subject_x, adni_calibration],
np.r_[person_sites, np.full(len(adni_calibration), 2)])
subject_x = combat.transform(subject_x, person_sites)
X_train = combat.transform(X_train, site_train)
X_test = combat.transform(X_test, site_test)
def transform_cohort(fc, site):
return scaler.transform(combat.transform(fc, np.full(len(fc), site)))
scaler = StandardScaler().fit(subject_x)
subject_x = scaler.transform(subject_x)
X_train = scaler.transform(X_train)
t_p_vector = np.zeros((4950, 2))
for i in range(4950):
t_p_vector[i, 0], t_p_vector[i, 1] = t_test(subject_x[subject_y==1, i], subject_x[subject_y==0, i], 0.05, False)
t_p_vector_all.append(t_p_vector)
X_train = X_train[:,t_p_vector[:,1]<=p_thre]
# Preset common parameters; no tuning against external results.
# dx arrays are float; bincount and indexing need integer labels.
sample_weights = (1.0 / np.bincount(inverse)[inverse]
/ np.bincount(subject_y.astype(int), minlength=2)[y_train.astype(int)])
sample_weights *= len(sample_weights) / sample_weights.sum()
if X_train.shape[1] == 0:
raise ValueError('No edges passed the training-only filter')
clf = LinearSVC(penalty='l1', C=0.1, max_iter=20000, dual=False,
random_state=9).fit(X_train, y_train, sample_weight=sample_weights)
pred_tr = np.expand_dims(clf.predict(X_train), -1)
weight = np.zeros((len(t_p_vector)))
weight[t_p_vector[:,1]<=p_thre] = clf.coef_.squeeze()
# weight[t_p_vector[:,1]<=p_thre] = clf.feature_importances_.squeeze()
weights.append(weight.squeeze())
X_test = scaler.transform(X_test)
half_fc_adni_correct_ = transform_cohort(half_fc_adni_correct, 2)
half_fc_adni_correct_logic_define_ = transform_cohort(sa_dem_fc_adni_new_define, 2)
half_fc_adni_other_correct_ = transform_cohort(half_fc_adni_other_correct, 2)
half_fc_adni_other_correct_ = half_fc_adni_other_correct_[:,t_p_vector[:,1]<=p_thre]
pred_adni_other = clf.decision_function(half_fc_adni_other_correct_)
other_adni_pred_label.append(pred_adni_other)
half_fc_oas_other_correct_ = transform_cohort(half_fc_oas_other_correct, 0)
half_fc_oas_other_correct_ = half_fc_oas_other_correct_[:,t_p_vector[:,1]<=p_thre]
pred_oas_other = clf.decision_function(half_fc_oas_other_correct_)
other_oas_pred_label.append(pred_oas_other)
half_fc_habs_other_correct_ = transform_cohort(half_fc_habs_other_correct, 20)
half_fc_habs_other_correct_ = half_fc_habs_other_correct_[:,t_p_vector[:,1]<=p_thre]
pred_habs_other = clf.decision_function(half_fc_habs_other_correct_)
other_habs_pred_label.append(pred_habs_other)
X_test = X_test[:,t_p_vector[:,1]<=p_thre]
predPROB_te = clf.decision_function(X_test)
pred_test_uni = []
predPROB_test_uni = []
for te in te_sub_unqiu:
predPROB = np.mean(predPROB_te[te_subs == te],0)
predPROB_test_uni.append(predPROB)
if predPROB>0:
pred_test_uni.append(1)
else:
pred_test_uni.append(0)
pred_test_uni = np.array(pred_test_uni)
predPROB_test_uni = np.array(predPROB_test_uni)
run_pred = (predPROB_te > 0).astype(int)
sen, spe = sens_spec(run_pred, y_test)
perform.append([accuracy_score(y_test, run_pred), sen, spe,
roc_auc_score(y_test, predPROB_te)])
subject_y, subject_margin = participant_scores(te_subs, y_test, predPROB_te)
subject_pred = (subject_margin > 0).astype(int)
sen, spe = sens_spec(subject_pred, subject_y)
perform_sub.append([accuracy_score(subject_y, subject_pred), sen, spe,
roc_auc_score(subject_y, subject_margin)])
sen, spe = sens_spec(pred_test_uni[site_test_te==0], y_test[te_id][site_test_te==0])
acc = accuracy_score(y_test[te_id][site_test_te==0],pred_test_uni[site_test_te==0])
AUC = roc_auc_score(y_test[te_id][site_test_te==0], predPROB_test_uni[site_test_te==0])
print('training acc oas: {}'.format(accuracy_score(y_train[site_train==0],pred_tr[site_train==0])))
perform_oas.append([acc, sen, spe, AUC])
acc = accuracy_score(y_test[te_id][site_test_te==20],pred_test_uni[site_test_te==20])
print('training acc hasb: {}'.format(accuracy_score(y_train[site_train==20],pred_tr[site_train==20])))
perform_habs.append(acc)
#################independent test
pred_ind = np.expand_dims(clf.predict(half_fc_adni_correct_[:,t_p_vector[:,1]<=p_thre]), -1)
predPROB_ind = clf.decision_function(half_fc_adni_correct_[:,t_p_vector[:,1]<=p_thre])
indep_pred_prob.append(predPROB_ind)
sen, spe = sens_spec(pred_ind, sa_dem_dx_adni)
acc = accuracy_score(sa_dem_dx_adni,pred_ind)
AUC = roc_auc_score(sa_dem_dx_adni, np.array(predPROB_ind))
indep_performance_all.append([acc, sen, spe, AUC])
#################independent test logic mem
pred_ind = np.expand_dims(clf.predict(half_fc_adni_correct_logic_define_[:,t_p_vector[:,1]<=p_thre]), -1)
predPROB_ind = clf.decision_function(half_fc_adni_correct_logic_define_[:,t_p_vector[:,1]<=p_thre])
indep_pred_prob_logic.append(predPROB_ind)
sen, spe = sens_spec(pred_ind, sa_dem_dx_adni_new_define)
acc = accuracy_score(sa_dem_dx_adni_new_define,pred_ind)
AUC = roc_auc_score(sa_dem_dx_adni_new_define, np.array(predPROB_ind))
indep_performance_logic.append([acc, sen, spe, AUC])
#################independent test cn and dem
half_fc_cn_dem_pair_ = transform_cohort(half_fc_cn_dem_pair, 2)
half_fc_cn_sa_pair_ = transform_cohort(half_fc_cn_sa_pair, 2)
predPROB_ind = clf.decision_function(half_fc_cn_dem_pair_[:,t_p_vector[:,1]<=p_thre])
indep_cn_dem_prob.append(predPROB_ind)
predPROB_ind = clf.decision_function(half_fc_cn_sa_pair_[:,t_p_vector[:,1]<=p_thre])
indep_cn_sa_prob.append(predPROB_ind)
other_adni_pred_label = np.array(other_adni_pred_label)
other_oas_pred_label = np.array(other_oas_pred_label)
other_habs_pred_label = np.array(other_habs_pred_label)
#################averaging runs
other_oas_pred_label_uni = []
half_fc_oas_other_correct_uni = []
idx_oas_uni = []
for i, sub in enumerate(oasis_sa_sub_sess_uniq):
other_oas_pred_label_uni.append(other_oas_pred_label[:,oasis_sa_sub_sess_uniq == sub].mean(-1))
half_fc_oas_other_correct_uni.append(half_fc_oas_other_correct[oasis_sa_sub_sess_uniq == sub,:].mean(0))
idx_oas_uni.append(i)
other_oas_pred_label_uni = np.array(other_oas_pred_label_uni)
half_fc_oas_other_correct_uni = np.array(half_fc_oas_other_correct_uni)
idx_oas_uni = np.array(idx_oas_uni)
perform_oas = np.array(perform_oas)
print(perform_oas.mean(0))
perform_habs = np.array(perform_habs)
print(perform_habs.mean(0))
perform = np.array(perform)
print("Discovery input-row level [acc, sen, spe, auc]:", perform.mean(0))
perform_sub = np.array(perform_sub)
print("Discovery participant level [acc, sen, spe, auc]:", perform_sub.mean(0))
indep_pred_prob = np.array(indep_pred_prob)
indep_pred_prob_mean = np.mean(indep_pred_prob, 0).squeeze()
indep_pred_all_mean = np.zeros((len(indep_pred_prob_mean)))
indep_pred_all_mean[indep_pred_prob_mean>0] = 1
indep_pred_all_mean[indep_pred_prob_mean<0] = 0
acc = accuracy_score(sa_dem_dx_adni,indep_pred_all_mean)
AUC = roc_auc_score(sa_dem_dx_adni, np.array(indep_pred_prob_mean))
sen, spe = sens_spec(indep_pred_all_mean, sa_dem_dx_adni)
adni_perf = [acc, sen, spe, AUC]
print('ensemble input-row level acc:{} sen:{} spe:{} auc:{}'.format(acc, sen, spe, AUC))
_, sa_dem_sub_adni_unqi = np.unique(sa_dem_sub_adni, return_index=True)
adni_y_sub, adni_margin_sub = participant_scores(
sa_dem_sub_adni, sa_dem_dx_adni, indep_pred_prob_mean)
adni_pred_sub = (adni_margin_sub > 0).astype(int)
acc = accuracy_score(adni_y_sub, adni_pred_sub)
AUC = roc_auc_score(adni_y_sub, adni_margin_sub)
sen, spe = sens_spec(adni_pred_sub, adni_y_sub)
adni_perf_sublevel = [acc, sen, spe, AUC]
print('ADNI ensemble participant level acc:{} sen:{} spe:{} auc:{}'.format(acc, sen, spe, AUC))
adni_perf_uni = [acc, sen, spe, AUC]
indep_pred_prob = np.array(indep_cn_dem_prob)
indep_pred_prob_mean = np.mean(indep_pred_prob, 0).squeeze()
indep_pred_all_mean = np.zeros((len(indep_pred_prob_mean)))
indep_pred_all_mean[indep_pred_prob_mean>0] = 1
indep_pred_all_mean[indep_pred_prob_mean<0] = 0
acc = accuracy_score(dx_cn_dem_pair,indep_pred_all_mean)
AUC = roc_auc_score(dx_cn_dem_pair, np.array(indep_pred_prob_mean))
sen, spe = sens_spec(indep_pred_all_mean, dx_cn_dem_pair)
adni_perf = [acc, sen, spe, AUC]
print('ensemble acc cn vs dem:{} sen:{} spe:{} auc:{}'.format(acc, sen, spe, AUC))
indep_pred_prob = np.array(indep_cn_sa_prob)
indep_pred_prob_mean = np.mean(indep_pred_prob, 0).squeeze()
indep_pred_all_mean = np.zeros((len(indep_pred_prob_mean)))
indep_pred_all_mean[indep_pred_prob_mean>0] = 1
indep_pred_all_mean[indep_pred_prob_mean<0] = 0
acc = accuracy_score(dx_cn_sa_pair,indep_pred_all_mean)
AUC = roc_auc_score(dx_cn_sa_pair, np.array(indep_pred_prob_mean))
sen, spe = sens_spec(indep_pred_all_mean, dx_cn_sa_pair)
adni_perf = [acc, sen, spe, AUC]
print('ensemble acc cn vs sa:{} sen:{} spe:{} auc:{}'.format(acc, sen, spe, AUC))
# # #################association of PCD in ADNI
other_sub_adni_uni, other_idx_adni_uni = np.unique(other_sub_adni, return_index=True)
other_sub_adni_unique = []
other_sess_adni_unique = []
other_idx_adni_unique = []
for sub in other_sub_adni_uni:
# sub = '002_S_20.50 '
mask = other_sub_adni == sub
sess_sub = other_sess_adni[mask]
idx = np.where(mask)[0]
if sum(mask) == 1:
other_sub_adni_unique.append(sub)
other_sess_adni_unique.append(other_sess_adni[mask][0])
other_idx_adni_unique.append(idx[0])
else:
sess_num = []
for i,sess in enumerate(sess_sub):
if 'bl' in sess:
other_sub_adni_unique.append(sub)
other_sess_adni_unique.append(other_sess_adni[mask][i])
other_idx_adni_unique.append(idx[i])
skip = True
break
else:
sess_num.append(float(sess.split('m')[-1]))
skip = False
if skip == False:
sess_num = np.array(sess_num)
other_sub_adni_unique.append(sub)
other_sess_adni_unique.append(other_sess_adni[mask][np.argmin(sess_num)])
other_idx_adni_unique.append(idx[np.argmin(sess_num)])
skip = True
other_sub_adni_unique = np.array(other_sub_adni_unique)
other_sess_adni_unique = np.array(other_sess_adni_unique)
other_idx_adni_unique = np.array(other_idx_adni_unique)
PCD_adni_interest = {}
adni_path = r'E:\PHD\learning\research\AD\data\ADNI'
adni_demo = pd.read_csv(os.path.join(adni_path, 'ADNIMERGE_01May2024.csv'))
adni_demo_sub = adni_demo['PTID']
adni_demo_sess = adni_demo['VISCODE']
adni_demo_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_demo_sub == sub.strip())&(adni_demo_sess == sess.split(' ')[0])) > 0:
adni_demo_interest.append(adni_demo[['PTEDUCAT', 'PTRACCAT', 'APOE4', 'FDG', 'PIB', 'AV45', 'FBB',
'ABETA', 'TAU', 'PTAU', 'CDRSB', 'ADAS11', 'ADAS13', 'MMSE', 'RAVLT_immediate',
'RAVLT_learning', 'RAVLT_forgetting', 'RAVLT_perc_forgetting', 'LDELTOTAL', 'DIGITSCOR',
'TRABSCOR', 'FAQ', 'MOCA', 'EcogPtMem', 'EcogPtLang', 'EcogPtVisspat', 'EcogPtPlan',
'EcogPtOrgan', 'EcogPtDivatt', 'EcogPtTotal', 'EcogSPMem', 'EcogSPLang', 'EcogSPVisspat',
'EcogSPPlan', 'EcogSPOrgan', 'EcogSPDivatt', 'EcogSPTotal']][
(adni_demo_sub == sub.strip())&(adni_demo_sess == sess.split(' ')[0])].values[0,:])
else:
adni_demo_interest.append([np.nan for i in range(len(['PTEDUCAT', 'PTRACCAT', 'APOE4', 'FDG', 'PIB', 'AV45', 'FBB',
'ABETA', 'TAU', 'PTAU', 'CDRSB', 'ADAS11', 'ADAS13', 'MMSE', 'RAVLT_immediate',
'RAVLT_learning', 'RAVLT_forgetting', 'RAVLT_perc_forgetting', 'LDELTOTAL', 'DIGITSCOR',
'TRABSCOR', 'FAQ', 'MOCA', 'EcogPtMem', 'EcogPtLang', 'EcogPtVisspat', 'EcogPtPlan',
'EcogPtOrgan', 'EcogPtDivatt', 'EcogPtTotal', 'EcogSPMem', 'EcogSPLang', 'EcogSPVisspat',
'EcogSPPlan', 'EcogSPOrgan', 'EcogSPDivatt', 'EcogSPTotal']))])
adni_demo_interest = np.array(adni_demo_interest)
adni_demo_interest[adni_demo_interest[:,1]!='White',1] = 0
adni_demo_interest[adni_demo_interest[:,1]=='White',1] = 1
adni_demo_interest[adni_demo_interest[:,7]=='>1700',7] = 2000
adni_demo_interest[adni_demo_interest[:,9]=='<8',9] = 8
adni_demo_interest[adni_demo_interest[:,17]<-300,17] = np.nan
PCD_adni_interest.update({"AGE": other_age_adni[other_idx_adni_unique], 'GENDER': other_sex_adni[other_idx_adni_unique], 'dx': other_dx_adni[other_idx_adni_unique],
'PTEDUCAT': adni_demo_interest[:,0],
'PTRACCAT': adni_demo_interest[:,1], 'APOE4': adni_demo_interest[:,2], 'FDG': adni_demo_interest[:,3],
'AV45': adni_demo_interest[:,5], 'FBB': adni_demo_interest[:,6], 'ABETA': adni_demo_interest[:,7],
'TAU': adni_demo_interest[:,8], 'PTAU': adni_demo_interest[:,9], 'CDRSB': adni_demo_interest[:,10],
'ADAS11': adni_demo_interest[:,11], 'ADAS13': adni_demo_interest[:,12], 'MMSE': adni_demo_interest[:,13],
'RAVLT_immediate': adni_demo_interest[:,14], 'RAVLT_learning': adni_demo_interest[:,15], 'RAVLT_forgetting': adni_demo_interest[:,16],
'RAVLT_perc_forgetting': adni_demo_interest[:,17], 'LDELTOTAL': adni_demo_interest[:,18],
'TRABSCOR': adni_demo_interest[:,20], 'FAQ': adni_demo_interest[:,21], 'MOCA': adni_demo_interest[:,22],
'EcogPtMem': adni_demo_interest[:,23], 'EcogPtLang': adni_demo_interest[:,24], 'EcogPtVisspat': adni_demo_interest[:,25],
'EcogPtPlan': adni_demo_interest[:,26], 'EcogPtOrgan': adni_demo_interest[:,27], 'EcogPtDivatt': adni_demo_interest[:,28],
'EcogPtTotal': adni_demo_interest[:,29], 'EcogSPMem': adni_demo_interest[:,30], 'EcogSPLang': adni_demo_interest[:,31],
'EcogSPVisspat': adni_demo_interest[:,32], 'EcogSPPlan': adni_demo_interest[:,33], 'EcogSPOrgan': adni_demo_interest[:,34],
'EcogSPDivatt': adni_demo_interest[:,35], 'EcogSPTotal': adni_demo_interest[:,36]})
adni_cdr = pd.read_csv(os.path.join(adni_path, 'CDR_07May2024.csv'))
adni_cdr_sub = adni_cdr['PTID']
adni_cdr_sess = adni_cdr['VISCODE2']
adni_cdr_sess[adni_cdr_sess == 'sc'] = 'bl'
adni_cdr_sess[adni_cdr_sess == 'f'] = 'bl'
adni_cdr_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_cdr_sub == sub.strip())&(adni_cdr_sess == sess.split(' ')[0])) > 0:
adni_cdr_interest.append(adni_cdr[['CDMEMORY', 'CDORIENT', 'CDJUDGE', 'CDCOMMUN', 'CDHOME', 'CDCARE', 'CDGLOBAL']][
(adni_cdr_sub == sub.strip())&(adni_cdr_sess == sess.split(' ')[0])].values[0,:])
else:
adni_cdr_interest.append([np.nan for i in range(len(['CDMEMORY', 'CDORIENT', 'CDJUDGE', 'CDCOMMUN', 'CDHOME', 'CDCARE', 'CDGLOBAL']))])
adni_cdr_interest = np.array(adni_cdr_interest)
PCD_adni_interest.update({'CDMEMORY': adni_cdr_interest[:,0], 'CDORIENT': adni_cdr_interest[:,1], 'CDJUDGE': adni_cdr_interest[:,2],
'CDCOMMUN': adni_cdr_interest[:,3], 'CDHOME': adni_cdr_interest[:,4], 'CDCARE': adni_cdr_interest[:,5],
'CDGLOBAL': adni_cdr_interest[:,6]})
adni_faq = pd.read_csv(os.path.join(adni_path, 'FAQ_08May2024.csv'))
adni_faq_sub = adni_faq['PTID']
adni_faq_sess = adni_faq['VISCODE2']
adni_faq_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_faq_sub == sub.strip())&(adni_faq_sess == sess.split(' ')[0])) > 0:
adni_faq_interest.append(adni_faq[['FAQFINAN', 'FAQFORM', 'FAQSHOP', 'FAQGAME', 'FAQBEVG', 'FAQMEAL',
'FAQEVENT', 'FAQTV', 'FAQREM', 'FAQTRAVL']][
(adni_faq_sub == sub.strip())&(adni_faq_sess == sess.split(' ')[0])].values[0,:])
else:
adni_faq_interest.append([np.nan for i in range(len(['FAQFINAN', 'FAQFORM', 'FAQSHOP', 'FAQGAME', 'FAQBEVG', 'FAQMEAL',
'FAQEVENT', 'FAQTV', 'FAQREM', 'FAQTRAVL']))])
adni_faq_interest = np.array(adni_faq_interest)
adni_faq_interest[(adni_faq_interest==1)|(adni_faq_interest==2)|(adni_faq_interest==3)] = 1
adni_faq_interest[(adni_faq_interest==4)] = 2
adni_faq_interest[(adni_faq_interest==5)] = 3
PCD_adni_interest.update({'FAQFINAN': adni_faq_interest[:,0], 'FAQFORM': adni_faq_interest[:,1], 'FAQSHOP': adni_faq_interest[:,2],
'FAQGAME': adni_faq_interest[:,3], 'FAQBEVG': adni_faq_interest[:,4], 'FAQMEAL': adni_faq_interest[:,5],
'FAQEVENT': adni_faq_interest[:,6], 'FAQTV': adni_faq_interest[:,7], 'FAQREM': adni_faq_interest[:,8],
'FAQTRAVL': adni_faq_interest[:,9]})
adni_neuroBIO = pd.read_csv(os.path.join(adni_path, 'UPENNBIOMK_ROCHE_ELECSYS_16May2024.csv'))
adni_neuroBIO_sub = adni_neuroBIO['RID']
adni_neuroBIO_sess = adni_neuroBIO['VISCODE2']
adni_neuroBIO_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_neuroBIO_sub == int(sub.strip().split("_")[-1]))&(adni_neuroBIO_sess == sess.split(' ')[0])) > 0:
adni_neuroBIO_interest.append(adni_neuroBIO[['ABETA42', 'TAU', 'PTAU']][
(adni_neuroBIO_sub == int(sub.strip().split("_")[-1]))&(adni_neuroBIO_sess == sess.split(' ')[0])].values[0,:])
else:
adni_neuroBIO_interest.append([np.nan for i in range(len(['ABETA42', 'TAU', 'PTAU']))])
adni_neuroBIO_interest = np.array(adni_neuroBIO_interest, dtype=object)
# These columns replace the ADNIMERGE ones cleaned above, so apply the same censoring convention here.
adni_neuroBIO_interest[adni_neuroBIO_interest[:,0]=='>1700',0] = 2000
adni_neuroBIO_interest[adni_neuroBIO_interest[:,2]=='<8',2] = 8
try:
adni_neuroBIO_interest = adni_neuroBIO_interest.astype(float)
except ValueError as exc:
raise ValueError('Unhandled censored CSF value in the Elecsys file; add it to the convention above') from exc
PCD_adni_interest.update({'ABETA': adni_neuroBIO_interest[:,0], 'TAU': adni_neuroBIO_interest[:,1], 'PTAU': adni_neuroBIO_interest[:,2]})
adni_gds = pd.read_csv(os.path.join(adni_path, 'GDSCALE_08May2024.csv'))
adni_gds_sub = adni_gds['PTID']
adni_gds_sess = adni_gds['VISCODE2']
adni_gds_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_gds_sub == sub.strip())&(adni_gds_sess == sess.split(' ')[0])) > 0:
adni_gds_interest.append(adni_gds['GDTOTAL'][
(adni_gds_sub == sub.strip())&(adni_gds_sess == sess.split(' ')[0])].values[0])
else:
adni_gds_interest.append(np.nan)
adni_gds_interest = np.array(adni_gds_interest)
PCD_adni_interest.update({'GDTOTAL': adni_gds_interest})
adni_his = pd.read_csv(os.path.join(adni_path, 'MODHACH_08May2024.csv'))
adni_his_sub = adni_his['PTID']
adni_his_sess = adni_his['VISCODE2']
adni_his_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_his_sub == sub.strip())) > 0:
adni_his_interest.append(adni_his[['HMSTEPWS', 'HMSOMATC', 'HMEMOTIO', 'HMHYPERT', 'HMSTROKE', 'HMNEURSM',
'HMNEURSG', 'HMSCORE', 'HMONSET']][
(adni_his_sub == sub.strip())].values[0,:])
else:
adni_his_interest.append([np.nan for i in range(len(['HMSTEPWS', 'HMSOMATC', 'HMEMOTIO', 'HMHYPERT', 'HMSTROKE', 'HMNEURSM',
'HMNEURSG', 'HMSCORE', 'HMONSET']))])
adni_his_interest = np.array(adni_his_interest)
PCD_adni_interest.update({'HMSTEPWS': adni_his_interest[:,0], 'HMSOMATC': adni_his_interest[:,1], 'HMEMOTIO': adni_his_interest[:,2],
'HMHYPERT': adni_his_interest[:,3], 'HMSTROKE': adni_his_interest[:,4], 'HMNEURSM': adni_his_interest[:,5],
'HMNEURSG': adni_his_interest[:,6], 'HMSCORE': adni_his_interest[:,7], 'HMONSET': adni_his_interest[:,8]})
adni_moca = pd.read_csv(os.path.join(adni_path, 'MOCA_08May2024.csv'))
adni_moca_sub = adni_moca['PTID']
adni_moca_sess = adni_moca['VISCODE2']
adni_moca_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_moca_sub == sub.strip())&(adni_moca_sess == sess.split(' ')[0])) > 0:
adni_moca_interest.append(adni_moca[['TRAILS', 'CUBE', 'CLOCKCON', 'CLOCKNO', 'CLOCKHAN', 'LION', 'RHINO', 'CAMEL',
'IMMT1W1', 'IMMT1W2', 'IMMT1W3', 'IMMT1W4', 'IMMT1W5', 'IMMT2W1', 'IMMT2W2', 'IMMT2W3',
'IMMT2W4', 'IMMT2W5', 'DIGFOR', 'DIGBACK', 'LETTERS', 'SERIAL1', 'SERIAL2', 'SERIAL3', 'SERIAL4',
'SERIAL5', 'REPEAT1', 'REPEAT2', 'FFLUENCY', 'ABSTRAN', 'ABSMEAS', 'DELW1', 'DELW2', 'DELW3',
'DELW4', 'DELW5', 'DATE', 'MONTH', 'YEAR', 'DAY', 'PLACE', 'CITY']][
(adni_moca_sub == sub.strip())&(adni_moca_sess == sess.split(' ')[0])].values[0,:])
else:
adni_moca_interest.append([np.nan for i in range(len(['TRAILS', 'CUBE', 'CLOCKCON', 'CLOCKNO', 'CLOCKHAN', 'LION', 'RHINO', 'CAMEL',
'IMMT1W1', 'IMMT1W2', 'IMMT1W3', 'IMMT1W4', 'IMMT1W5', 'IMMT2W1', 'IMMT2W2', 'IMMT2W3',
'IMMT2W4', 'IMMT2W5', 'DIGFOR', 'DIGBACK', 'LETTERS', 'SERIAL1', 'SERIAL2', 'SERIAL3', 'SERIAL4',
'SERIAL5', 'REPEAT1', 'REPEAT2', 'FFLUENCY', 'ABSTRAN', 'ABSMEAS', 'DELW1', 'DELW2', 'DELW3',
'DELW4', 'DELW5', 'DATE', 'MONTH', 'YEAR', 'DAY', 'PLACE', 'CITY']))])
adni_moca_interest = np.array(adni_moca_interest)
adni_moca_subtract = adni_moca_interest[:,21:26].sum(-1)
adni_moca_subtract[(adni_moca_subtract>=2) & (adni_moca_subtract<=3)] = 2
adni_moca_subtract[adni_moca_subtract>=4] = 3
adni_moca_frq = adni_moca_interest[:,28]
adni_moca_frq[adni_moca_frq<11] = 0
adni_moca_frq[adni_moca_frq>=11] = 1
PCD_adni_interest.update({'visuospatial_executive': adni_moca_interest[:,:5].sum(-1), 'naming': adni_moca_interest[:,5:8].sum(-1),
'abstraction': adni_moca_interest[:,29:31].sum(-1), 'delayed_recall': adni_moca_interest[:,31:36].sum(-1)/3,
'orientation': adni_moca_interest[:,36:].sum(-1), 'attention': adni_moca_interest[:,18:20].sum(-1) + adni_moca_interest[:,20] + adni_moca_subtract,
'language': adni_moca_interest[:,26:28].sum(-1) + adni_moca_frq})
adni_nps = pd.read_csv(os.path.join(adni_path, 'NPIQ_08May2024.csv'))
adni_nps_sub = adni_nps['PTID']
adni_nps_sess = adni_nps['VISCODE2']
adni_nps_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_nps_sub == sub.strip())&(adni_nps_sess == sess.split(' ')[0])) > 0:
adni_nps_interest.append(adni_nps[['NPIASEV', 'NPIBSEV', 'NPICSEV', 'NPIDSEV', 'NPIESEV', 'NPIFSEV', 'NPIGSEV',
'NPIHSEV', 'NPIISEV', 'NPIJSEV', 'NPIKSEV', 'NPILSEV','NPISCORE']][
(adni_nps_sub == sub.strip())&(adni_nps_sess == sess.split(' ')[0])].values[0,:])
else:
adni_nps_interest.append([np.nan for i in range(len(['NPIASEV', 'NPIBSEV', 'NPICSEV', 'NPIDSEV', 'NPIESEV', 'NPIFSEV', 'NPIGSEV',
'NPIHSEV', 'NPIISEV', 'NPIJSEV', 'NPIKSEV', 'NPILSEV','NPISCORE']))])
adni_nps_interest = np.array(adni_nps_interest)
adni_nps_interest[adni_nps_interest == -4] = 0
adni_nps_interest[np.isnan(adni_nps_interest)] = 0
PCD_adni_interest.update({'NPIASEV': adni_nps_interest[:,0], 'NPIBSEV': adni_nps_interest[:,1], 'NPICSEV': adni_nps_interest[:,2],
'NPIDSEV': adni_nps_interest[:,3], 'NPIESEV': adni_nps_interest[:,4], 'NPIFSEV': adni_nps_interest[:,5],
'NPIGSEV': adni_nps_interest[:,6], 'NPIHSEV': adni_nps_interest[:,7], 'NPIISEV': adni_nps_interest[:,8],
'NPIJSEV': adni_nps_interest[:,9], 'NPIKSEV': adni_nps_interest[:,10], 'NPILSEV': adni_nps_interest[:,11],
'NPISCORE': adni_nps_interest[:,12]})
adni_neu = pd.read_csv(os.path.join(adni_path, 'NEUROBAT_08May2024.csv'))
adni_neu_sub = adni_neu['PTID']
adni_neu_sess = adni_neu['VISCODE2']
adni_neu_interest = []
for sub, sess in zip(other_sub_adni_unique, other_sess_adni_unique):
if sum((adni_neu_sub == sub.strip())&(adni_neu_sess == sess.split(' ')[0])) > 0:
adni_neu_interest.append(adni_neu[['LIMMTOTAL', 'DSPANFOR', 'DSPANFLTH', 'DSPANBAC', 'DSPANBLTH', 'CATANIMSC', 'CATANPERS',
'CATVEGESC', 'CATVGPERS', 'TRAASCOR', 'TRABSCOR', 'DIGITSCOR','LDELTOTAL', 'BNTTOTAL']][
(adni_neu_sub == sub.strip())&(adni_neu_sess == sess.split(' ')[0])].values[0,:])
else:
adni_neu_interest.append([np.nan for i in range(len(['LIMMTOTAL', 'DSPANFOR', 'DSPANFLTH', 'DSPANBAC', 'DSPANBLTH', 'CATANIMSC',
'CATANPERS', 'CATVEGESC', 'CATVGPERS', 'TRAASCOR', 'TRABSCOR', 'DIGITSCOR','LDELTOTAL',
'BNTTOTAL']))])
adni_neu_interest = np.array(adni_neu_interest)
adni_neu_interest[adni_neu_interest == -4] = 0
PCD_adni_interest.update({'LIMMTOTAL': adni_neu_interest[:,0], 'CATANIMSC': adni_neu_interest[:,5],
'CATANPERS': adni_neu_interest[:,6], 'CATVEGESC': adni_neu_interest[:,7], 'DSPANFOR': adni_neu_interest[:,1],
'TRAASCOR': adni_neu_interest[:,9], 'DSPANBAC': adni_neu_interest[:,3], 'DIGITSCOR': adni_neu_interest[:,11],
'BNTTOTAL': adni_neu_interest[:,13]})
# ---- Phenotype group comparisons ----
PCD_keys = ['AGE', 'GENDER', 'PTEDUCAT', 'PTRACCAT', 'APOE4', 'FDG', 'AV45', 'FBB', 'ABETA', 'TAU', 'PTAU', 'CDRSB',
'ADAS11', 'ADAS13', 'MMSE', 'RAVLT_immediate', 'RAVLT_learning', 'RAVLT_forgetting',
'RAVLT_perc_forgetting', 'LDELTOTAL', 'TRABSCOR', 'FAQ', 'MOCA', 'CDMEMORY',
'CDORIENT', 'CDJUDGE', 'CDCOMMUN', 'CDHOME', 'CDCARE', 'CDGLOBAL', 'FAQFINAN', 'FAQFORM',
'FAQSHOP', 'FAQGAME', 'FAQBEVG', 'FAQMEAL', 'FAQEVENT', 'FAQTV', 'FAQREM', 'FAQTRAVL', 'GDTOTAL',
'visuospatial_executive', 'naming', 'abstraction', 'delayed_recall', 'orientation', 'attention',
'language', 'NPIASEV', 'NPIBSEV', 'NPICSEV',
'NPIDSEV', 'NPIESEV', 'NPIFSEV', 'NPIGSEV', 'NPIHSEV', 'NPIISEV', 'NPIJSEV', 'NPIKSEV', 'NPILSEV',
'NPISCORE', 'LIMMTOTAL', 'DSPANFOR', 'DSPANBAC', 'CATANIMSC', 'CATANPERS', 'CATVEGESC', 'TRAASCOR', 'DIGITSCOR', 'BNTTOTAL']
df_adni = pd.DataFrame(PCD_adni_interest).astype(float)
df_adni['A+'] = (df_adni['ABETA'] <= float(os.environ['ARC_CSF_ABETA_CUTOFF'])) * 1.0
df_adni.loc[df_adni['ABETA'].isnull(), 'A+'] = np.nan
df_adni['T+'] = (df_adni['PTAU'] > float(os.environ['ARC_CSF_PTAU_CUTOFF'])) * 1.0
df_adni.loc[df_adni['PTAU'].isnull(), 'T+'] = np.nan
PCD_keys.append('A+')
PCD_keys.append('T+')
stat_all_adni_cn = []
stat_keys_adni = []
stat_test = []
mask_cn = (other_dx_adni[other_idx_adni_unique]==2)
mask_MCI = (other_dx_adni[other_idx_adni_unique]==3)
# other_adni_pred_label_avg = other_adni_pred_label.mean(0)[other_idx_adni_unique]
other_adni_pred_label_mask = other_adni_pred_label>0
other_adni_pred_label_avg = other_adni_pred_label_mask.mean(0)[other_idx_adni_unique]
for i, key in enumerate(PCD_keys):
stat_keys_adni.append(key)
fea = df_adni[key]
fea_demlike_cn = fea[mask_cn&(other_adni_pred_label_avg>=0.5)]
fea_demlike_cn = fea_demlike_cn[np.isfinite(fea_demlike_cn)]
fea_salike_cn = fea[mask_cn&(other_adni_pred_label_avg<0.5)]
fea_salike_cn = fea_salike_cn[np.isfinite(fea_salike_cn)]
fea_unstable_cn = fea[mask_cn&(other_adni_pred_label_avg!=0)&(other_adni_pred_label_avg!=1)]
fea_unstable_cn = fea_unstable_cn[np.isfinite(fea_unstable_cn)]
if len(fea_demlike_cn)>5 and len(fea_salike_cn)>5 and (len(fea_salike_cn) + len(fea_demlike_cn))>20 and (len(set(fea_salike_cn))>1 or len(set(fea_demlike_cn)) > 1):
if key == 'A+' or key == 'T+':
[F_cn,p_cn,_,_]=chi2_contingency([[np.sum(fea_salike_cn==1), np.sum(fea_demlike_cn==1)],
[np.sum(fea_salike_cn==0), np.sum(fea_demlike_cn==0)]])
cohens_d = cohen_d(fea_salike_cn, fea_demlike_cn)
stat_all_adni_cn.append([F_cn, p_cn, cohens_d, len(fea_demlike_cn), len(fea_salike_cn)])
if key == 'GENDER' or key == 'PTRACCAT':
plt.figure(figsize=(5,10))
ax = plt.gca()
count = np.sum(fea_salike_cn==1)
ax.bar(0, count, color = red, width = 1, fill=False, edgecolor= red, linewidth=5, hatch = '+')
count = np.sum(fea_demlike_cn==1)
ax.bar(1, count, color = blue, width = 1, fill=False, edgecolor= blue, linewidth=5, hatch = '+')
if key == 'PTRACCAT':
count = np.sum(fea_salike_cn==0)
else:
count = np.sum(fea_salike_cn==2)
ax.bar(2, count, color = red, width = 1, fill=False, edgecolor= red, linewidth=5, hatch = 'x')
if key == 'PTRACCAT':
count = np.sum(fea_demlike_cn==0)
else:
count = np.sum(fea_demlike_cn==2)
ax.bar(3, count, color = blue, width = 1, fill=False, edgecolor= blue, linewidth=5, hatch = 'x')
ax.set_yticks((0, 100, 200, 300, 400))
ax.set_xticks([1.5])
ax.set_xticklabels(['CN'])
ax.set_xlabel(key, **font2)
ax.spines['top'].set_color('none') # 设置上‘脊梁’为红色
ax.spines['right'].set_color('none') # 设置上‘脊梁’为无色
else:
plt.figure(figsize=(5,10))
ax = plt.gca()
count = np.sum(fea_demlike_cn<0.5)
ax.bar(0, count, color = red, width = 1, fill=False, edgecolor= red, linewidth=5, hatch = '+')
count = np.sum(fea_salike_cn<0.5)
ax.bar(1, count, color = blue, width = 1, fill=False, edgecolor= blue, linewidth=5, hatch = '+')
count = np.sum(fea_demlike_cn==1)
ax.bar(2, count, color = red, width = 1, fill=False, edgecolor= red, linewidth=5, hatch = 'x')
count = np.sum(fea_salike_cn==1)
ax.bar(3, count, color = blue, width = 1, fill=False, edgecolor= blue, linewidth=5, hatch = 'x')
ax.set_yticks((0, 50, 100, 150))
ax.set_xticks([1.5])
ax.set_xticklabels(['CN'])
ax.set_xlabel(key, **font2)
ax.spines['top'].set_color('none') # 设置上‘脊梁’为红色
ax.spines['right'].set_color('none') # 设置上‘脊梁’为无色
else:
F_cn,p_cn=kruskal(fea_salike_cn, fea_demlike_cn)
cohens_d = cohen_d(fea_salike_cn, fea_demlike_cn)
stat_all_adni_cn.append([F_cn, p_cn, cohens_d, len(fea_demlike_cn), len(fea_salike_cn)])
result1 = np.r_[fea_salike_cn, fea_demlike_cn]
dx_label = np.r_[['CN' for i in range(len(fea_salike_cn)+len(fea_demlike_cn))]]
cluster_label = np.r_[['SA-like' for i in range(len(fea_salike_cn))], ['DEM-like' for i in range(len(fea_demlike_cn))]]
df = pd.DataFrame({'data': result1, 'DX': dx_label, 'Cluster': cluster_label})
plt.figure(figsize=(5,10))
ax = plt.gca()
my_pal = {'SA-like': blue, 'DEM-like': red}
# my_pal = {'SA-like': blue, 1: red, 2: 'limegreen', 3: RGB(0,.8,0.1)}
sns.stripplot(data=df, x="DX", y="data", hue="Cluster", palette = my_pal, ax = ax, dodge=True, legend = False, size = 5, jitter=0.2)
sns.violinplot(data=df, x="DX", y="data", hue="Cluster", palette = my_pal, fill=False, ax = ax, saturation=1, inner = 'box')
ax.scatter([-0.2,0.2],y= np.r_[fea_salike_cn.median(), fea_demlike_cn.median()],c="white",s = 100, edgecolor = 'grey', zorder=3, linewidth = 2)
plt.xlim(-1, 1)
for collection in ax.collections:
if isinstance(collection, matplotlib.collections.PolyCollection):
collection.set_edgecolor(collection.get_facecolor())
collection.set_facecolor('none')
for h in ax.legend_.legendHandles:
if isinstance(h, matplotlib.patches.Rectangle):
h.set_edgecolor(h.get_facecolor())
h.set_facecolor('none')
h.set_linewidth(2.5)
ax.get_legend().remove()
unit = (np.max(fea[np.isfinite(fea)]) - np.min(fea[np.isfinite(fea)]))//2
if abs(unit) < 1:
unit = round((np.max(fea[np.isfinite(fea)]) - np.min(fea[np.isfinite(fea)]))/2, 2)
ax.set_yticks(np.arange(np.min(fea[np.isfinite(fea)])-unit//2, np.max(fea[np.isfinite(fea)])+unit//2, unit))
else:
ax.set_yticks(np.arange(np.min(fea[np.isfinite(fea)])-unit//2, np.max(fea[np.isfinite(fea)])+unit//2, unit).astype(int))
if key == 'AGE':
ax.set_yticks(np.arange(20, 110, 30))
plt.ylim(20, 110)
ax.set_ylabel(key)
ax.spines['top'].set_color('none') # 设置上‘脊梁’为红色
ax.spines['right'].set_color('none') # 设置上‘脊梁’为无色
# plt.savefig(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\cluster\stats_compare\0.5\{}_adni_f.svg'.format(key), bbox_inches = 'tight')
else:
stat_all_adni_cn.append([np.nan, np.nan, np.nan, len(fea_demlike_cn), len(fea_salike_cn)])
stat_all_adni_cn = np.array(stat_all_adni_cn)
stat_all_adni_cn[np.r_[[70,71], np.arange(5,11)],1]
out = multipletests(stat_all_adni_cn[np.r_[[70,71], np.arange(5,11)],1], method="fdr_bh")[1]
np.array(stat_keys_adni)[np.r_[[70,71], np.arange(5,11)]]
stat_adni_nonan_fdrp = multipletests(stat_all_adni_cn[np.isfinite(stat_all_adni_cn[:,1]),1], method="fdr_bh")[1]
stat_keys_adni_nonan_fdrp = np.array(stat_keys_adni)[np.isfinite(stat_all_adni_cn[:,1])]
stat_adni_nonan = stat_all_adni_cn[np.isfinite(stat_all_adni_cn[:,1])]
# sio.savemat(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\cluster\stats_compare\0.5\stats_compare_ini_adni.mat', {'stat_keys_adni':stat_keys_adni,
# 'stat_all_adni_cn': stat_all_adni_cn})
# #################association of PCD in OASIS
oasis_sa_sub = pd.Series(oasis_sa_sub_sess_uniq).str.split('_', expand = True).iloc[:,0]
oasis_sa_sub_unique, unique_id = np.unique(oasis_sa_sub, return_index=True)
oasis_sa_sub_sess_unique = oasis_sa_sub_sess_uniq[unique_id]
PCD_all_interest = {}
oasis_demo = pd.read_csv(os.path.join(oasis_path, 'UDS-A2-InformatDemos.csv'))
oasis_demo2 = pd.read_csv(os.path.join(oasis_path, 'OASIS3_UDSa2_cs_demo.csv'))
oasis_demo_sub = oasis_demo['Subject']
oasis_demo_sub2 = oasis_demo2['OASISID']
sa_dem_sex_oas = []
sa_dem_age_oas = []
for i, sub in enumerate(sa_dem_sub_oas_):
# sa_dem_sex_oas.append(sex_oas[(sub_oas==sub)&(sess_oas == sa_dem_sess_oas[i])])
if sum(oasis_demo_sub == sub) >= 1:
sa_dem_sex_oas.append(oasis_demo['INSEX'][oasis_demo_sub == sub].values[0])
else:
if sum(oasis_demo_sub2 == sub) >= 1:
sa_dem_sex_oas.append(oasis_demo2['INSEX'][oasis_demo_sub2 == sub].values[0])
sa_dem_age_oas.append(age_oas[(sub_oas==sub)&(sess_oas == sa_dem_sess_oas[i])])
sa_dem_sex_oas = np.array(sa_dem_sex_oas).squeeze() # female: 2; male: 1
sa_dem_age_oas = np.array(sa_dem_age_oas).squeeze()
oasis_path = r'E:\PHD\learning\research\AD\data\OASIS3'
oasis_demo = pd.read_csv(os.path.join(oasis_path, 'UDS-A2-InformatDemos.csv'))
oasis_demo2 = pd.read_csv(os.path.join(oasis_path, 'OASIS3_UDSa2_cs_demo.csv'))
oasis_demo_sub = oasis_demo['Subject']
oasis_demo_sub2 = oasis_demo2['OASISID']
oasis_demo_interest = []
for sub in oasis_sa_sub_sess_unique:
sub = sub[:8]
if sum(oasis_demo_sub == sub) >= 1:
oasis_demo_interest.append(oasis_demo[['INRACE', 'INEDUC', 'INSEX']][oasis_demo_sub == sub].values[0,:])
else:
if sum(oasis_demo_sub2 == sub) >= 1:
oasis_demo_interest.append(oasis_demo2[['INRACE', 'INEDUC', 'INSEX']][oasis_demo_sub2 == sub].values[0,:])
else:
oasis_demo_interest.append([np.nan, np.nan, np.nan])
oasis_demo_interest = np.array(oasis_demo_interest)
PCD_all_interest.update({"INRACE": oasis_demo_interest[:,0], 'INEDUC': oasis_demo_interest[:,1], 'age': other_age_oas[unique_id],
'sex': oasis_demo_interest[:,2], 'dx': other_dx_oas_encode[unique_id]})
oas_pcd_info = sio.loadmat(r'E:\PHD\learning\research\AD\data\OASIS3\fMRIdata_OASIS3_struct.mat')
tvar = oas_pcd_info['OASIS3_Phenotypic_struct']
PCD_data = []
pcd_keys = ['subjectID_Date', 'apoe', 'NPIQINF', 'NPIQINFX', 'DEL',
'DELSEV', 'HALL', 'HALLSEV', 'AGIT', 'AGITSEV', 'DEPD', 'DEPDSEV', 'ANX',
'ANXSEV', 'ELAT', 'ELATSEV', 'APA', 'APASEV', 'DISN', 'DISNSEV', 'IRR',
'IRRSEV', 'MOT', 'MOTSEV', 'NITE', 'NITESEV', 'APP', 'APPSEV', 'mmse',
'cdr', 'commun','homehobb', 'judgment', 'memory', 'orient', 'perscare', 'sumbox',
'DIGIF', 'DIGIB', 'ANIMALS', 'VEG', 'TRAILA', 'TRAILALI', 'TRAILB',
'TRAILBLI', 'WAIS', 'LOGIMEM', 'MEMUNITS', 'MEMTIME', 'BOSTON']
for idx in range(len(tvar)):
sub_data = []
for name in pcd_keys:
value = tvar[idx][0][name]
for i in range(3):
if isinstance(value, np.ndarray):
try:
value = value[0]
except IndexError:
if name == 'DX':
value = 'nan'
else:
value = "error"
else:
if name == 'subjectID_Date' or name == 'dx1':
value = str(value)
sub_data.append(value)
break
PCD_data.append(sub_data)
PCD_data = np.array(PCD_data)
PCD_data[:,1] = pd.Series(PCD_data[:,1]).str.count('4')
#update NPI 0 score
for i in [5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27]:
mask = PCD_data[:, i] == 'nan'
PCD_data[np.isfinite(PCD_data[:,2].astype(float))&mask, i] = 0
PCD_data[np.isfinite(PCD_data[:,2].astype(float))][:,[5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27]] = 0
PCD_sub = np.array([sub[:8] for sub in PCD_data[:,0]])
PCD_day = np.array([sub[10:] for sub in PCD_data[:,0]]).astype(int)
PCD_sess = []
day_min = []
for da in PCD_day:
mask = day_all==da
if sum(mask) > 0:
PCD_sess.append(sess_all[day_all==da][0])
else:
diff = day_all-da
diff_min_id = np.argmin(abs(diff))
diff_min = diff[diff_min_id]
day_min.append(diff_min)
PCD_sess.append(sess_all[diff_min_id])
PCD_sess = np.array(PCD_sess)
PCD_interest = []
for sub_sess in oasis_sa_sub_sess_unique:
sub = sub_sess[:8]
sess = sub_sess[9:]
if sum((PCD_sub == sub) & (PCD_sess == sess))== 1:
PCD_interest.append(PCD_data[:,np.r_[[1, 2, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25], np.arange(27,50)]][
(PCD_sub == sub) & (PCD_sess == sess)].squeeze())
elif sum((PCD_sub == sub) & (PCD_sess == sess))> 1:
pcd = PCD_data[:,np.r_[[1, 2, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25], np.arange(27,50)]][
(PCD_sub == sub) & (PCD_sess == sess)].squeeze()
pcd_new = []
for i in range(pcd.shape[-1]):
info = pcd[:,i]
if np.isnan(info.astype(float)).sum() == len(info):
pcd_new.append(np.nan)
else:
pcd_new.append(info[~np.isnan(info.astype(float))][0].squeeze())
PCD_interest.append(np.array(pcd_new))
else:
PCD_interest.append(np.array([np.nan for i in range(len(['apoe', 'NPIQINF', 'DELSEV', 'HALLSEV', 'AGITSEV', 'DEPDSEV',
'ANXSEV', 'ELATSEV', 'APASEV', 'DISNSEV', 'IRRSEV', 'MOTSEV', 'NITESEV',
'APPSEV', 'mmse', 'cdr', 'commun','homehobb', 'judgment', 'memory', 'orient', 'perscare', 'sumbox',
'DIGIF', 'DIGIB', 'ANIMALS', 'VEG', 'TRAILA', 'TRAILALI', 'TRAILB',
'TRAILBLI', 'WAIS', 'LOGIMEM', 'MEMUNITS', 'MEMTIME', 'BOSTON']))]))
PCD_interest = np.array(PCD_interest)
PCD_all_interest.update({"apoe": PCD_interest[:,0], 'NPIQINF': PCD_interest[:,1], 'DELSEV': PCD_interest[:,2],
'HALLSEV': PCD_interest[:,3], 'AGITSEV': PCD_interest[:,4], 'DEPDSEV': PCD_interest[:,5],
'ANXSEV': PCD_interest[:,6], 'ELATSEV': PCD_interest[:,7], 'APASEV': PCD_interest[:,8],
'DISNSEV': PCD_interest[:,9], 'IRRSEV': PCD_interest[:,10], 'MOTSEV': PCD_interest[:,11],
'NITESEV': PCD_interest[:,12], 'APPSEV': PCD_interest[:,13], 'mmse': PCD_interest[:,14],
'cdr': PCD_interest[:,15], 'commun': PCD_interest[:,16], 'homehobb': PCD_interest[:,17],
'judgment': PCD_interest[:,18], 'memory': PCD_interest[:,19], 'orient': PCD_interest[:,20], 'perscare': PCD_interest[:,21],
'sumbox': PCD_interest[:,22], 'DIGIF': PCD_interest[:,23], 'DIGIB': PCD_interest[:,24],
'ANIMALS': PCD_interest[:,25], 'VEG': PCD_interest[:,26], 'TRAILA': PCD_interest[:,27],
'TRAILALI': PCD_interest[:,28], 'TRAILB': PCD_interest[:,29], 'TRAILBLI': PCD_interest[:,30], 'WAIS': PCD_interest[:,31],
'LOGIMEM': PCD_interest[:,32], 'MEMUNITS': PCD_interest[:,33], 'MEMTIME': PCD_interest[:,34], 'BOSTON': PCD_interest[:,35]})
oasis_faq = pd.read_csv(os.path.join(oasis_path, 'UDS-B7-FAQ.csv'))
oasis_faq_sub_day = oasis_faq['UDS_B7FAQDATA ID']
oasis_faq_sub = oasis_faq_sub_day.str.split('_', expand = True).iloc[:,0]
oasis_faq_day = oasis_faq_sub_day.str.split('_d', expand = True).iloc[:,-1].astype(int)
oasis_faq_sess = []
day_min = []
for da in oasis_faq_day:
mask = day_all==da
if sum(mask) > 0:
oasis_faq_sess.append(sess_all[day_all==da][0])
else:
diff = day_all-da
diff_min_id = np.argmin(abs(diff))
diff_min = diff[diff_min_id]
day_min.append(diff_min)
oasis_faq_sess.append(sess_all[diff_min_id])
oasis_faq_sess = np.array(oasis_faq_sess)
oasis_faq_interest = []
for sub_sess in oasis_sa_sub_sess_unique:
sub = sub_sess[:8]
sess = sub_sess[9:]
if sum((oasis_faq_sub == sub) & (oasis_faq_sess == sess))> 0:
oasis_faq_interest.append(oasis_faq[['BILLS', 'TAXES', 'SHOPPING', 'GAMES', 'STOVE', 'MEALPREP', 'EVENTS',
'PAYATTN', 'REMDATES', 'TRAVEL']][
(oasis_faq_sub == sub) & (oasis_faq_sess == sess)].values[0,:])
else:
oasis_faq_interest.append([np.nan for i in range(len(['BILLS', 'TAXES', 'SHOPPING', 'GAMES', 'STOVE', 'MEALPREP', 'EVENTS',
'PAYATTN', 'REMDATES', 'TRAVEL']))])
oasis_faq_interest = np.array(oasis_faq_interest)
oasis_faq_interest[(oasis_faq_interest==8)|(oasis_faq_interest==9)] = np.nan
PCD_all_interest.update({"BILLS": oasis_faq_interest[:,0], 'TAXES': oasis_faq_interest[:,1], 'SHOPPING': oasis_faq_interest[:,2],
'GAMES': oasis_faq_interest[:,3], 'STOVE': oasis_faq_interest[:,4], 'MEALPREP': oasis_faq_interest[:,5],
'EVENTS': oasis_faq_interest[:,6], 'PAYATTN': oasis_faq_interest[:,7], 'REMDATES': oasis_faq_interest[:,8],
'TRAVEL': oasis_faq_interest[:,9]})
oasis_pet = pd.read_csv(os.path.join(oasis_path, 'OASIS3_PUP.csv'))
oasis_pet['Centil_fBP_TOT_CORTMEAN_av45'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fBP_TOT_CORTMEAN_av45'] = oasis_pet['Centil_fBP_TOT_CORTMEAN'][oasis_pet['tracer'] == 'AV45']
oasis_pet['Centil_fSUVR_TOT_CORTMEAN_av45'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fSUVR_TOT_CORTMEAN_av45'] = oasis_pet['Centil_fSUVR_TOT_CORTMEAN'][oasis_pet['tracer'] == 'AV45']
oasis_pet['Centil_fBP_rsf_TOT_CORTMEAN_av45'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fBP_rsf_TOT_CORTMEAN_av45'] = oasis_pet['Centil_fBP_rsf_TOT_CORTMEAN'][oasis_pet['tracer'] == 'AV45']
oasis_pet['Centil_fSUVR_rsf_TOT_CORTMEAN_av45'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fSUVR_rsf_TOT_CORTMEAN_av45'] = oasis_pet['Centil_fSUVR_rsf_TOT_CORTMEAN'][oasis_pet['tracer'] == 'AV45']
oasis_pet['Centil_fBP_TOT_CORTMEAN_pib'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fBP_TOT_CORTMEAN_pib'] = oasis_pet['Centil_fBP_TOT_CORTMEAN'][oasis_pet['tracer'] == 'PIB']
oasis_pet['Centil_fSUVR_TOT_CORTMEAN_pib'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fSUVR_TOT_CORTMEAN_pib'] = oasis_pet['Centil_fSUVR_TOT_CORTMEAN'][oasis_pet['tracer'] == 'PIB']
oasis_pet['Centil_fBP_rsf_TOT_CORTMEAN_pib'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fBP_rsf_TOT_CORTMEAN_pib'] = oasis_pet['Centil_fBP_rsf_TOT_CORTMEAN'][oasis_pet['tracer'] == 'PIB']
oasis_pet['Centil_fSUVR_rsf_TOT_CORTMEAN_pib'] = np.zeros((oasis_pet.shape[0]))
oasis_pet['Centil_fSUVR_rsf_TOT_CORTMEAN_pib'] = oasis_pet['Centil_fSUVR_rsf_TOT_CORTMEAN'][oasis_pet['tracer'] == 'PIB']
oasis_pet_sub_day = oasis_pet['PUP_PUPTIMECOURSEDATA ID']
oasis_pet_sub = oasis_pet_sub_day.str.split('_', expand = True).iloc[:,0]
oasis_pet_day = oasis_pet_sub_day.str.split('_d', expand = True).iloc[:,-1].astype(int)
oasis_pet_sess = []
day_min = []
for da in oasis_pet_day:
mask = day_all==da
if sum(mask) > 0:
oasis_pet_sess.append(sess_all[day_all==da][0])
else:
diff = day_all-da
diff_min_id = np.argmin(abs(diff))
diff_min = diff[diff_min_id]
day_min.append(diff_min)
oasis_pet_sess.append(sess_all[diff_min_id])
oasis_pet_sess = np.array(oasis_pet_sess)
oasis_pet_interest = []
for sub_sess in oasis_sa_sub_sess_unique:
sub = sub_sess[:8]
sess = sub_sess[9:]
if sum((oasis_pet_sub == sub) & (oasis_pet_sess == sess))> 0:
oasis_pet_interest.append(oasis_pet[['Centil_fBP_TOT_CORTMEAN_av45', 'Centil_fSUVR_TOT_CORTMEAN_av45',
'Centil_fBP_rsf_TOT_CORTMEAN_av45', 'Centil_fSUVR_rsf_TOT_CORTMEAN_av45',
'Centil_fBP_TOT_CORTMEAN_pib', 'Centil_fSUVR_TOT_CORTMEAN_pib',
'Centil_fBP_rsf_TOT_CORTMEAN_pib', 'Centil_fSUVR_rsf_TOT_CORTMEAN_pib']][
(oasis_pet_sub == sub) & (oasis_pet_sess == sess)].values[0,:])
else:
oasis_pet_interest.append([np.nan for i in range(len(['Centil_fBP_TOT_CORTMEAN_av45', 'Centil_fSUVR_TOT_CORTMEAN_av45',
'Centil_fBP_rsf_TOT_CORTMEAN_av45', 'Centil_fSUVR_rsf_TOT_CORTMEAN_av45',
'Centil_fBP_TOT_CORTMEAN_pib', 'Centil_fSUVR_TOT_CORTMEAN_pib',
'Centil_fBP_rsf_TOT_CORTMEAN_pib', 'Centil_fSUVR_rsf_TOT_CORTMEAN_pib']))])
oasis_pet_interest = np.array(oasis_pet_interest)
PCD_all_interest.update({"Centil_fBP_TOT_CORTMEAN_av45": oasis_pet_interest[:,0], 'Centil_fSUVR_TOT_CORTMEAN_av45': oasis_pet_interest[:,1],
'Centil_fBP_rsf_TOT_CORTMEAN_av45': oasis_pet_interest[:,2], 'Centil_fSUVR_rsf_TOT_CORTMEAN_av45': oasis_pet_interest[:,3],
"Centil_fBP_TOT_CORTMEAN_pib": oasis_pet_interest[:,4], 'Centil_fSUVR_TOT_CORTMEAN_pib': oasis_pet_interest[:,5],
'Centil_fBP_rsf_TOT_CORTMEAN_pib': oasis_pet_interest[:,6], 'Centil_fSUVR_rsf_TOT_CORTMEAN_pib': oasis_pet_interest[:,7]})
PCD_all_interest['TRAILB'] = PCD_all_interest['TRAILB'].astype(float)
PCD_all_interest['TRAILB'][PCD_all_interest['TRAILB']>400] = np.nan
df_oas = pd.DataFrame(PCD_all_interest).astype(float)
PCD_keys = ['INRACE', 'INEDUC', 'sex', 'age', 'apoe', 'DELSEV', 'HALLSEV', 'AGITSEV',
'DEPDSEV', 'ANXSEV', 'ELATSEV', 'APASEV', 'DISNSEV', 'IRRSEV', 'MOTSEV', 'NITESEV',
'APPSEV', 'mmse', 'cdr', 'commun','homehobb', 'judgment', 'memory', 'orient', 'perscare', 'sumbox', 'VEG',
'DIGIF', 'DIGIB', 'ANIMALS', 'TRAILA', 'TRAILALI', 'TRAILB', 'TRAILBLI', 'WAIS', 'LOGIMEM',
'MEMUNITS', 'MEMTIME', 'BOSTON', 'BILLS', 'TAXES', 'SHOPPING', 'GAMES', 'STOVE', 'MEALPREP',
'EVENTS', 'PAYATTN', 'REMDATES', 'TRAVEL', 'Centil_fBP_TOT_CORTMEAN_av45', 'Centil_fSUVR_TOT_CORTMEAN_av45',
'Centil_fBP_rsf_TOT_CORTMEAN_av45', 'Centil_fSUVR_rsf_TOT_CORTMEAN_av45', 'Centil_fBP_TOT_CORTMEAN_pib', 'Centil_fSUVR_TOT_CORTMEAN_pib',
'Centil_fBP_rsf_TOT_CORTMEAN_pib', 'Centil_fSUVR_rsf_TOT_CORTMEAN_pib']
stat_all_oas_cn = []
stat_keys_oas = []
stat_test = []
# mask_cn = (other_dx_oas_encode==2)&(other_age_oas>=75)
mask_MCI = (other_dx_oas_encode[unique_id]==3)
# mask_cn = ((other_dx_oas_encode[unique_id]==2)|(other_dx_oas_encode[unique_id]==3))
mask_cn = (other_dx_oas_encode[unique_id]==2)
other_oas_pred_label_mask = other_oas_pred_label>0
other_oas_pred_label_avg = other_oas_pred_label_mask.mean(0)[unique_id]
# other_oas_pred_label_avg = other_oas_pred_label.mean(0)[unique_id]
oasis_sa_sub_sess_uniq_ = pd.Series(oasis_sa_sub_sess_unique).str.split('_', expand = True).iloc[:,0]
# #################association of PCD in HABS sa
half_fc_habs_aging_correct_uniq = []
sub_ind_ = pd.Series(sub_ind).str.split('_', expand = True).iloc[:,0]
other_sub_habs = pd.Series(other_sub_sess_habs).str.split('_', expand = True).iloc[:,0]
habs_sub_uni, uni_idx = np.unique(other_sub_habs, return_index=True)
age_habs_aging_correct_uniq = []
sex_habs_aging_correct_uniq = []
for i in range(len(habs_sub_uni)):
half_fc_habs_aging_correct_uniq.append(half_fc_habs_other_correct[other_sub_habs == habs_sub_uni[i], :].mean(0))
age_habs_aging_correct_uniq.append(other_age_habs[other_sub_habs == habs_sub_uni[i]].mean(0))
sex_habs_aging_correct_uniq.append(other_sex_habs[other_sub_habs == habs_sub_uni[i]][0])
half_fc_habs_aging_correct_uniq = np.array(half_fc_habs_aging_correct_uniq)
age_habs_aging_correct_uniq = np.array(age_habs_aging_correct_uniq)
sex_habs_aging_correct_uniq = np.array(sex_habs_aging_correct_uniq)
habs_aging_behavior_path = r'I:\data\habs'
cog_path = os.path.join(habs_aging_behavior_path, 'Cognition_HABS_DataRelease_2.0.csv')
cog_info = pd.read_csv(cog_path)
demo_path = os.path.join(habs_aging_behavior_path, 'Demographics_HABS_DataRelease_2.0.csv')
demo_info = pd.read_csv(demo_path)
demo_info = demo_info.rename(columns={'SubjID': 'SubjIDshort'})
clinic_path = os.path.join(habs_aging_behavior_path, 'ClinicalMeasures_HABS_DataRelease_2.0.csv')
clinic_info = pd.read_csv(clinic_path)
FDG_path = os.path.join(habs_aging_behavior_path, 'FDG_FS6_SUVR_HABS_DataRelease_2.0.csv')
FDG_info = pd.read_csv(FDG_path)
FDG_info['fdg_mean'] = FDG_info.iloc[:,4:].mean(1)
FTP_path = os.path.join(habs_aging_behavior_path, 'FTP_FS6_SUVR_HABS_DataRelease_2.0.csv')
FTP_info = pd.read_csv(FTP_path)
FTP_info['ftp_mean'] = FTP_info.iloc[:,4:].mean(1)
PIB_path = os.path.join(habs_aging_behavior_path, 'PIB_FS6_DVR_HABS_DataRelease_2.0.csv')
PIB_info = pd.read_csv(PIB_path)
PIB_info.loc[PIB_info['StudyArc'] == 'HAB_18m_PET', 'StudyArc'] = 'HAB_3.0'
FDG_info.loc[FDG_info['StudyArc'] == 'HAB_18m_PET', 'StudyArc'] = 'HAB_3.0'
FDG_info.loc[FDG_info['StudyArc'] == 'HAB_AI', 'StudyArc'] = 'HAB_6.0'
pcd_all_habs = pd.merge(cog_info, demo_info, on=['SubjIDshort', 'StudyArc', 'NP_SessionDate'], how='outer')
pcd_all_habs = pd.merge(pcd_all_habs, clinic_info, on=['SubjIDshort', 'StudyArc', 'NP_SessionDate'], how='outer')
pcd_all_habs = pd.merge(pcd_all_habs, FDG_info, on=['SubjIDshort', 'StudyArc'], how='outer')
pcd_all_habs = pd.merge(pcd_all_habs, FTP_info, on=['SubjIDshort', 'StudyArc'], how='outer')
pcd_all_habs = pd.merge(pcd_all_habs, PIB_info, on=['SubjIDshort', 'StudyArc'], how='outer')
pcd_all_habs['sub_sess'] = pcd_all_habs['SubjIDshort'].str.split('_', expand = True).iloc[:,-1] + '_' +pcd_all_habs['StudyArc']
# Build both masks before recoding; match OASIS INRACE (1 = White), non-White = 2, missing stays NaN.
race_raw = pcd_all_habs['Race'].copy()
pcd_all_habs['Race'] = np.where(race_raw.isnull(), np.nan, np.where(race_raw == 'W', 1.0, 2.0))
_, hpc_fc_idx, pcd_idx = np.intersect1d(habs_sub_uni, pd.Series(pcd_all_habs['sub_sess']).str.split('_', expand = True).iloc[:,0], return_indices=True)
PCD_habs_interest = {}
PCD_habs_interest['age']= other_age_habs[uni_idx][hpc_fc_idx]
PCD_habs_interest['sex']= other_sex_habs[uni_idx][hpc_fc_idx]
PCD_habs_interest['YrsOfEd']= pcd_all_habs['YrsOfEd'][pcd_idx]
PCD_habs_interest['APOE_haplotype']= pcd_all_habs['APOE_haplotype'][pcd_idx].astype(str).str.count('4')
PCD_habs_interest['BNT_30']= pcd_all_habs['BNT_30'][pcd_idx]
PCD_habs_interest['CAT3']= pcd_all_habs['CAT3'][pcd_idx]
PCD_habs_interest['CAT_Animal_Total']= pcd_all_habs['CAT_Animal_Total'][pcd_idx]
PCD_habs_interest['CAT_Vegetable_Total']= pcd_all_habs['CAT_Vegetable_Total'][pcd_idx]
PCD_habs_interest['CAT_Fruit_Total']= pcd_all_habs['CAT_Fruit_Total'][pcd_idx]
PCD_habs_interest['Digits_Forward']= pcd_all_habs['Digits_Forward'][pcd_idx]
PCD_habs_interest['Digits_Backwards']= pcd_all_habs['Digits_Backwards'][pcd_idx]
PCD_habs_interest['DigitSym']= pcd_all_habs['DigitSym'][pcd_idx]
PCD_habs_interest['FAS_Total']= pcd_all_habs['FAS_Total'][pcd_idx]
PCD_habs_interest['FAS_F_Total']= pcd_all_habs['FAS_F_Total'][pcd_idx]
PCD_habs_interest['FAS_A_Total']= pcd_all_habs['FAS_A_Total'][pcd_idx]
PCD_habs_interest['FAS_S_Total']= pcd_all_habs['FAS_S_Total'][pcd_idx]
PCD_habs_interest['FCsrt_FNC']= pcd_all_habs['FCsrt_FNC'][pcd_idx]
PCD_habs_interest['FCsrt_Free']= pcd_all_habs['FCsrt_Free'][pcd_idx]
PCD_habs_interest['LetterNum_Total']= pcd_all_habs['LetterNum_Total'][pcd_idx]
PCD_habs_interest['LogicMem_IL']= pcd_all_habs['LogicMem_IL'][pcd_idx]
PCD_habs_interest['LogicMem_DR']= pcd_all_habs['LogicMem_DR'][pcd_idx]
PCD_habs_interest['SRT_dr']= pcd_all_habs['SRT_dr'][pcd_idx]
PCD_habs_interest['SRT_cltr']= pcd_all_habs['SRT_cltr'][pcd_idx]
PCD_habs_interest['SRT_cr']= pcd_all_habs['SRT_cr'][pcd_idx]
PCD_habs_interest['SRT_ltr']= pcd_all_habs['SRT_ltr'][pcd_idx]
PCD_habs_interest['SRT_lts']= pcd_all_habs['SRT_lts'][pcd_idx]
PCD_habs_interest['SRT_mc']= pcd_all_habs['SRT_mc'][pcd_idx]
PCD_habs_interest['SRT_str']= pcd_all_habs['SRT_str'][pcd_idx]
PCD_habs_interest['SRT_tr']= pcd_all_habs['SRT_tr'][pcd_idx]
PCD_habs_interest['TMT_A']= pcd_all_habs['TMT_A'][pcd_idx]
PCD_habs_interest['TMT_B']= pcd_all_habs['TMT_B'][pcd_idx]
PCD_habs_interest['VFDT']= pcd_all_habs['VFDT'][pcd_idx]
PCD_habs_interest['CDR_Global']= pcd_all_habs['CDR_Global'][pcd_idx]
PCD_habs_interest['CDR_SB']= pcd_all_habs['CDR_SB'][pcd_idx]
PCD_habs_interest['CDR_Memory']= pcd_all_habs['CDR_Memory'][pcd_idx]
PCD_habs_interest['MMSE_Total']= pcd_all_habs['MMSE_Total'][pcd_idx]
PCD_habs_interest['MMSE_Orientation']= pcd_all_habs['MMSE_Orientation'][pcd_idx]
PCD_habs_interest['MMSE_ImmRecall']= pcd_all_habs['MMSE_ImmRecall'][pcd_idx]
PCD_habs_interest['MMSE_AttnCalc']= pcd_all_habs['MMSE_AttnCalc'][pcd_idx]
PCD_habs_interest['MMSE_DelRecall']= pcd_all_habs['MMSE_DelRecall'][pcd_idx]
PCD_habs_interest['MMSE_Language']= pcd_all_habs['MMSE_Language'][pcd_idx]
PCD_habs_interest['MMSE_Pentagons']= pcd_all_habs['MMSE_Pentagons'][pcd_idx]
PCD_habs_interest['GDS_Total']= pcd_all_habs['GDS_Total'][pcd_idx]
PCD_habs_interest['Hachinski']= pcd_all_habs['Hachinski'][pcd_idx]
PCD_habs_interest['Race']= pcd_all_habs['Race'][pcd_idx]
PCD_keys = ['Race', 'YrsOfEd', 'sex', 'age', 'APOE_haplotype', 'BNT_30', 'CAT3', 'CAT_Animal_Total', 'CAT_Vegetable_Total', 'CAT_Fruit_Total',
'Digits_Forward', 'Digits_Backwards', 'DigitSym', 'FAS_Total', 'FAS_F_Total', 'FAS_A_Total', 'FAS_S_Total', 'FCsrt_FNC', 'FCsrt_Free', 'LetterNum_Total',
'LogicMem_IL', 'LogicMem_DR', 'SRT_dr', 'SRT_cltr', 'SRT_cr', 'SRT_ltr', 'SRT_lts', 'SRT_mc', 'SRT_str', 'SRT_tr', 'TMT_A', 'TMT_B', 'VFDT', 'CDR_Global',
'CDR_SB', 'CDR_Memory', 'MMSE_Total', 'MMSE_Orientation', 'MMSE_ImmRecall', 'MMSE_AttnCalc', 'MMSE_DelRecall', 'MMSE_Language', 'MMSE_Pentagons',
'GDS_Total', 'Hachinski']
df_pcd_habs = pd.DataFrame(PCD_habs_interest).astype(float)
# ---- Phenotype group comparisons ----
stat_all_habs_cn = []
stat_keys_habs = []
stat_test = []
other_habs_pred_label_mask = other_habs_pred_label>0
other_habs_pred_label_avg = other_habs_pred_label_mask.mean(0)[uni_idx[hpc_fc_idx]]
df_all = {'age': np.r_[df_oas['age'], df_pcd_habs['age']],
# 'sub': np.r_[np.delete(HPC_sa_sub_uni,idx), oasis_sa_sub_sess_uniq_, habs_sub_uni[hpc_fc_idx]],
'sub': np.r_[oasis_sa_sub_sess_uniq_, habs_sub_uni[hpc_fc_idx]],
'sex': np.r_[df_oas['sex'], df_pcd_habs['sex']],
'cluster_label': np.r_[other_oas_pred_label_avg, other_habs_pred_label_avg],
# White = 1, non-White = 2 in both cohorts.
'race': np.r_[np.where(df_oas['INRACE'].isnull(), np.nan, np.where(df_oas['INRACE'] == 1, 1.0, 2.0)),
df_pcd_habs['Race']],
'educ': np.r_[df_oas['INEDUC'], df_pcd_habs['YrsOfEd']],
'dx': np.r_[df_oas['dx'], other_dx_habs[uni_idx[hpc_fc_idx]]],
'tmt_a': np.r_[df_oas['TRAILA'], df_pcd_habs['TMT_A']],
'tmt_b': np.r_[df_oas['TRAILB'], df_pcd_habs['TMT_B']],
# OASIS apoe is already an e4 count (see PCD_data[:,1] above).
'apoe': np.r_[df_oas['apoe'], df_pcd_habs['APOE_haplotype']],
'DELSEV': np.r_[df_oas['DELSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'HALLSEV': np.r_[df_oas['HALLSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'AGITSEV': np.r_[df_oas['AGITSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'DEPDSEV': np.r_[df_oas['DEPDSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'ANXSEV': np.r_[df_oas['ANXSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'ELATSEV': np.r_[df_oas['ELATSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'APASEV': np.r_[df_oas['APASEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'DISNSEV': np.r_[df_oas['DISNSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'IRRSEV': np.r_[df_oas['IRRSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'MOTSEV': np.r_[df_oas['MOTSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'NITESEV': np.r_[df_oas['NITESEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'APPSEV': np.r_[df_oas['APPSEV'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'mmse': np.r_[df_oas['mmse'], df_pcd_habs['MMSE_Total']],
'cdr': np.r_[df_oas['cdr'], df_pcd_habs['CDR_Global']],
'commun': np.r_[df_oas['commun'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'homehobb': np.r_[df_oas['homehobb'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'perscare': np.r_[df_oas['perscare'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'judgment': np.r_[df_oas['judgment'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'memory': np.r_[df_oas['memory'], df_pcd_habs['CDR_Memory']],
'orient': np.r_[df_oas['orient'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'sumbox': np.r_[df_oas['sumbox'], df_pcd_habs['CDR_SB']],
'RAVLT_immediate': np.r_[[np.nan for i in range(len(df_oas['sex']))], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'RAVLT_learning': np.r_[[np.nan for i in range(len(df_oas['sex']))], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'RAVLT_forgetting': np.r_[[np.nan for i in range(len(df_oas['sex']))], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'RAVLT_perc_forgetting': np.r_[[np.nan for i in range(len(df_oas['sex']))], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'LogMem_Delay': np.r_[df_oas['MEMUNITS'], df_pcd_habs['LogicMem_DR']],
'LogMem_immediate': np.r_[df_oas['LOGIMEM'], df_pcd_habs['LogicMem_IL']],
'CFT_ANIMALS': np.r_[df_oas['ANIMALS'], df_pcd_habs['CAT_Animal_Total']],
'CFT_VEG': np.r_[df_oas['VEG'], df_pcd_habs['CAT_Vegetable_Total']],
'DIGIF': np.r_[df_oas['DIGIF'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'DIGIB': np.r_[df_oas['DIGIB'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'BOSTON': np.r_[df_oas['BOSTON'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'BILLS': np.r_[df_oas['BILLS'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'TAXES': np.r_[df_oas['TAXES'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'SHOPPING': np.r_[df_oas['SHOPPING'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'GAMES': np.r_[df_oas['GAMES'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'STOVE': np.r_[df_oas['STOVE'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'MEALPREP': np.r_[df_oas['MEALPREP'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'EVENTS': np.r_[df_oas['EVENTS'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'PAYATTN': np.r_[df_oas['PAYATTN'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'REMDATES': np.r_[df_oas['REMDATES'], [np.nan for i in range(len(df_pcd_habs['sex']))]],
'TRAVEL': np.r_[df_oas['TRAVEL'], [np.nan for i in range(len(df_pcd_habs['sex']))]]}
df_all = pd.DataFrame(df_all)
stat_discovery_cn = []
stat_keys_discovery = []
stat_test = []
mask_cn = (df_all['dx']==2)#&(df_all['age']>=75)
PCD_all_keys = list(df_all.columns)
for i, key in enumerate(PCD_all_keys):
if key in {'sub', 'cluster_label', 'dx'}:
continue
stat_keys_discovery.append(key)
fea = df_all[key]
fea_demlike_cn = fea[mask_cn&(df_all['cluster_label']>=.5)]
fea_demlike_cn = fea_demlike_cn[np.isfinite(fea_demlike_cn)]
fea_salike_cn = fea[mask_cn&(df_all['cluster_label']<.5)]
fea_salike_cn = fea_salike_cn[np.isfinite(fea_salike_cn)]
fea_unstable_cn = fea[mask_cn&(df_all['cluster_label']!=0)&(df_all['cluster_label']!=1)]
fea_unstable_cn = fea_unstable_cn[np.isfinite(fea_unstable_cn)]
if len(fea_demlike_cn)>5 and len(fea_salike_cn)>5 and (len(fea_salike_cn) + len(fea_demlike_cn))>20 and (len(set(fea_salike_cn)) >1 or len(set(fea_demlike_cn)) > 1):
F_cn,p_cn=kruskal(fea_salike_cn, fea_demlike_cn)
if key == 'sex' or key == 'race':
[F_cn,p_cn,_,_]=chi2_contingency([[np.sum(fea_salike_cn==1), np.sum(fea_demlike_cn==1)],
[np.sum((fea_salike_cn==2)|(fea_salike_cn==0)), np.sum((fea_demlike_cn==2)|(fea_demlike_cn==0))]])
cohens_d = cohen_d(fea_salike_cn, fea_demlike_cn)
stat_discovery_cn.append([F_cn, p_cn, cohens_d, len(fea_demlike_cn), len(fea_salike_cn)])
if key == 'sex' or key == 'race':
plt.figure(figsize=(5,10))
ax = plt.gca()
count = np.sum(fea_salike_cn==1)
ax.bar(0, count, color = red, width = 1, fill=False, edgecolor= red, linewidth=5, hatch = '+')
count = np.sum(fea_demlike_cn==1)
ax.bar(1, count, color = blue, width = 1, fill=False, edgecolor= blue, linewidth=5, hatch = '+')
if key == 'race':
count = np.sum(fea_salike_cn==2)
else:
count = np.sum(fea_salike_cn==2)
ax.bar(2, count, color = red, width = 1, fill=False, edgecolor= red, linewidth=5, hatch = 'x')
if key == 'race':
count = np.sum(fea_demlike_cn==2)
else:
count = np.sum(fea_demlike_cn==2)
ax.bar(3, count, color = blue, width = 1, fill=False, edgecolor= blue, linewidth=5, hatch = 'x')
ax.set_yticks((0, 200, 400, 600, 800))
ax.set_xticks([1.5])
ax.set_xticklabels(['CN'])
ax.set_xlabel(key, **font2)
ax.spines['top'].set_color('none') # 设置上‘脊梁’为红色
ax.spines['right'].set_color('none') # 设置上‘脊梁’为无色
else:
result1 = np.r_[fea_salike_cn, fea_demlike_cn]
dx_label = np.r_[['CN' for i in range(len(fea_salike_cn)+len(fea_demlike_cn))]]
cluster_label = np.r_[['SA-like' for i in range(len(fea_salike_cn))], ['DEM-like' for i in range(len(fea_demlike_cn))]]
df = pd.DataFrame({'data': result1, 'DX': dx_label, 'Cluster': cluster_label})
plt.figure(figsize=(5,10))
ax = plt.gca()
my_pal = {'SA-like': blue, 'DEM-like': red}
# my_pal = {'SA-like': blue, 1: red, 2: 'limegreen', 3: RGB(0,.8,0.1)}
sns.stripplot(data=df, x="DX", y="data", hue="Cluster", palette = my_pal, ax = ax, dodge=True, legend = False, size = 5, jitter=0.2)
sns.violinplot(data=df, x="DX", y="data", hue="Cluster", palette = my_pal, fill=False, ax = ax, saturation=1, inner = 'box')
ax.scatter([-0.2,0.2],y= np.r_[fea_salike_cn.median(), fea_demlike_cn.median()],c="white",s = 100, edgecolor = 'grey', zorder=3, linewidth = 2)
plt.xlim(-1, 1)
for collection in ax.collections:
if isinstance(collection, matplotlib.collections.PolyCollection):
collection.set_edgecolor(collection.get_facecolor())
collection.set_facecolor('none')
for h in ax.legend_.legendHandles:
if isinstance(h, matplotlib.patches.Rectangle):
h.set_edgecolor(h.get_facecolor())
h.set_facecolor('none')
h.set_linewidth(2.5)
ax.get_legend().remove()
unit = (np.max(fea[np.isfinite(fea)]) - np.min(fea[np.isfinite(fea)]))//2
if abs(unit) < 1:
unit = round((np.max(fea[np.isfinite(fea)]) - np.min(fea[np.isfinite(fea)]))/2, 2)
ax.set_yticks(np.arange(np.min(fea[np.isfinite(fea)])-unit, np.max(fea[np.isfinite(fea)])+unit, unit))
else:
ax.set_yticks(np.arange(np.min(fea[np.isfinite(fea)])-unit, np.max(fea[np.isfinite(fea)])+unit, unit).astype(int))
if key == 'age':
ax.set_yticks(np.arange(20, 110, 30))
plt.ylim(20, 110)
ax.set_ylabel(key)
ax.spines['top'].set_color('none') # 设置上‘脊梁’为红色
ax.spines['right'].set_color('none') # 设置上‘脊梁’为无色
# plt.savefig(r'F:\PHD\learning\project\super_age\more_sub\SA_defined_ricado\cluster\stats_compare\0.5\{}_discovery_nodelete.svg'.format(key), bbox_inches = 'tight')
else:
stat_discovery_cn.append([np.nan, np.nan, np.nan, len(fea_demlike_cn), len(fea_salike_cn)])
stat_discovery_cn = np.array(stat_discovery_cn)
valid = np.isfinite(stat_discovery_cn[:, 1])
stat_discovery_cn_fdr = stat_discovery_cn.copy()
if valid.any():
stat_discovery_cn_fdr[valid, 1] = multipletests(stat_discovery_cn[valid, 1], method='fdr_bh')[1]
stat_keys_discovery_nonan = np.array(stat_keys_discovery)[valid]
###################################demograph in discovery
sub_ind_ = pd.Series(sub_ind).str.split('_', expand = True).iloc[:,0]
sub_ind_uni, sub_idx = np.unique(sub_ind_, return_index=True)
oasis_demo_interest_classification = []
for sub in sub_ind_uni:
if sum(oasis_demo_sub == sub) >= 1:
oasis_demo_interest_classification.append(oasis_demo[['INRACE', 'INEDUC']][oasis_demo_sub == sub].values[0,:])
else:
if sum(oasis_demo_sub2 == sub) >= 1:
oasis_demo_interest_classification.append(oasis_demo2[['INRACE', 'INEDUC']][oasis_demo_sub2 == sub].values[0,:])
else:
oasis_demo_interest_classification.append([np.nan, np.nan])
oasis_demo_interest_classification = np.array(oasis_demo_interest_classification)
habs_demo_interest_classification = []
habs_demo_sub = pcd_all_habs["SubjIDshort"].str.split('_', expand = True).iloc[:,-1]
for sub in sub_ind_uni:
if sum(habs_demo_sub == sub) >= 1:
habs_demo_interest_classification.append(pcd_all_habs[['Race', 'YrsOfEd']][habs_demo_sub == sub].values[0,:])
else:
habs_demo_interest_classification.append([np.nan, np.nan])
habs_demo_interest_classification = np.array(habs_demo_interest_classification)
oasis_demo_interest_classification[np.isnan(oasis_demo_interest_classification).sum(1)!=0] = habs_demo_interest_classification[np.isnan(oasis_demo_interest_classification).sum(1)!=0]
dx_ind_uni = dx_ind[sub_idx]
sex_ind_uni = sex_ind[sub_idx]
age_ind_uni = age_ind[sub_idx]
site_ind_uni = site_ind[sub_idx]
chi2_contingency([[np.sum(sex_ind_uni[(dx_ind_uni==0)]==1),np.sum(sex_ind_uni[(dx_ind_uni==1)]==1)],
[np.sum(sex_ind_uni[(dx_ind_uni==0)]==2),np.sum(sex_ind_uni[(dx_ind_uni==1)]==2)]])
# chi2_contingency([[np.sum(sex_ind_uni[(sex_ind_uni==0)&(site_ind_uni==1)]==1),np.sum(sex_ind_uni[(sex_ind_uni==1)&(site_ind_uni==1)]==1)],
# [np.sum(sex_ind_uni[(sex_ind_uni==0)&(site_ind_uni==1)]==2),np.sum(sex_ind_uni[(sex_ind_uni==1)&(site_ind_uni==1)]==2)]])
Q1 = np.percentile(age_ind_uni[dx_ind_uni==1], 25)
Q3 = np.percentile(age_ind_uni[dx_ind_uni==1], 75)
print(Q1)
print(Q3)
ranksums(age_ind_uni[dx_ind_uni==0], age_ind_uni[dx_ind_uni==1])
chi2_contingency([[np.sum(oasis_demo_interest_classification[(dx_ind_uni==0)]==1),np.sum(oasis_demo_interest_classification[(dx_ind_uni==1)]==1)],
[np.sum(oasis_demo_interest_classification[(dx_ind_uni==0)]==2),np.sum(oasis_demo_interest_classification[(dx_ind_uni==1)]==2)]])
# chi2_contingency([[np.sum(sex_ind_uni[(sex_ind_uni==0)&(site_ind_uni==1)]==1),np.sum(sex_ind_uni[(sex_ind_uni==1)&(site_ind_uni==1)]==1)],
# [np.sum(sex_ind_uni[(sex_ind_uni==0)&(site_ind_uni==1)]==2),np.sum(sex_ind_uni[(sex_ind_uni==1)&(site_ind_uni==1)]==2)]])
Q1 = np.percentile(oasis_demo_interest_classification[(dx_ind_uni==0)&np.isfinite(oasis_demo_interest_classification[:,1]),1], 25)
Q3 = np.percentile(oasis_demo_interest_classification[(dx_ind_uni==0)&np.isfinite(oasis_demo_interest_classification[:,1]),1], 75)
np.median(oasis_demo_interest_classification[(dx_ind_uni==0)&np.isfinite(oasis_demo_interest_classification[:,1]),1])
print(Q1)
print(Q3)
ranksums(oasis_demo_interest_classification[(dx_ind_uni==0)&np.isfinite(oasis_demo_interest_classification[:,1]),1], oasis_demo_interest_classification[(dx_ind_uni==1)&np.isfinite(oasis_demo_interest_classification[:,1]),1])
site_indx = []
for i,sub in enumerate(sub_ind_uni):
if 'HAB' in sub:
site_indx.append(i)
site_indx = np.array(site_indx)

clustering.py at commit ba44ae2, no license · at the source

Overview

Authors: Kanhao Zhao^1,2,3, Hua Xie⁴, Gregory A. Fonzo⁵, Nancy B. Carlisle⁶, Tovia Jacobs⁷, Ricardo S. Osorio⁸, Arpana Church³, Feng Vankee Lin^1,9,10, Yu Zhang^1,9,10, Alzheimer’s Disease Neuroimaging Initiative

ORCID iDs: Kanhao Zhao, Hua Xie, Gregory A. Fonzo, Nancy B. Carlisle, Arpana Church, Feng Vankee Lin, Yu Zhang

Department of Psychiatry and Behavioral Sciences, Stanford University School of Medicine, Stanford, CA, USA
Department of Bioengineering, Lehigh University, Bethlehem, PA, USA
Goodman-Luskin Microbiome Center, Vache and Tamar Manoukian Division of Digestive Diseases, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA, USA
Center for Neuroscience Research, Children’s National Hospital, Washington, DC, USA
Center for Psychedelic Research and Therapy, Department of Psychiatry and Behavioral Sciences, Dell Medical School, The University of Texas at Austin, Austin, TX, USA
Department of Psychology, Lehigh University, Bethlehem, PA, USA
Penn State College of Medicine, Hershey, PA, USA
Department of Psychiatry, New York University Grossman School of Medicine, New York, NY, USA
Wu Tsai Neurosciences Institute, Stanford University, Stanford, CA, USA
Stanford Institute for Human-Centered Artificial Intelligence, Stanford, CA, USA

Journal: Science advances, volume 12, issue 45, article eaed9351

Dates: received 15 November 2025; accepted 2 October 2026; published online 4 November 2026; in print November 2026

Type: Research article · Language: English

License: CC BY-NC

Identifiers: DOI · PMCID PMC13733855

Status: code verified (how the code is found and verified)

Categories: Alzheimer's / dementia (population)

Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging

Journal subjects: Neuroscience, Health and Medicine

Funding: National Science Foundation (2019035); national institutes of health (R21AG080425, R01MH129694, R01AG081768, R01MH132784); Brain & Behavior Research Foundation; Alzheimer's Association (AARG-22-972541); Stanford Knight Initiative for Brain Resilience and the Rosenkranz Foundation grant (RAR1-03); Lehigh University (CORE grant); One Mind - Baszucki Brain Research Fund; SEAL Future Fund

Citations: not cited yet (Europe PMC); 71 references in the paper

Abstract

Identifying the neurobiological basis of cognitive vulnerability is essential for preventing Alzheimer’s disease (AD). We derived the Alzheimer’s Resilience-Vulnerability Connectome (ARC) by contrasting SuperAgers (SA), an extreme phenotype of preserved late-life memory, with AD using resting-state fMRI. ARC distinguished SA from AD in a discovery cohort (N = 290; AUC = 0.85) and replicated independently (N = 143), involving prefrontal, temporal, and insular networks and associating with brain age. Applied to cognitively unimpaired adults (discovery: N = 818; replication: N = 497), ARC-based subtyping revealed SA-like and AD-like subgroups with comparable baseline cognition but divergent longitudinal trajectories. AD-like individuals exhibited faster cognitive decline, greater amyloid-β accumulation, and higher AD conversion risk. Exploratory GWAS nominated CLYBL and FRMD6 as candidate novel loci associated with the subtype difference. These findings suggest that ARC captures a continuous connectomic dimension of cognitive aging spanning SA, cognitively typical aging, and AD, highlighting the potential of functional connectomic variation for early risk stratification in AD.

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

Repositories

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

Zenodo 20073656

License: CC-BY-4.0

State: the link answers, verified on 7 November 2026

Evidence: files inventoried

Size: 0 files, 0 scripts

Software Heritage: not checked

Found in: “Data, code, and materials availability:”

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

Availability: 1 check, the latest on 7 November 2026: the link answers (HTTP 200)

7 November 2026: the link answers (HTTP 200)

At the source: zenodo.org/records/20073656

zhangyubrain/Alzheimers-resilient-connectome

License: none: the authors keep all their rights

State: the link answers, verified on 7 November 2026

Evidence: files inventoried

Commit: ba44ae26a27c55a13b3eb69fcdfae8361e73a5ba, 5 November 2026

Languages: Python (3)

Size: 4 files, 3 scripts

Software Heritage: not archived

Found in: “Data, code, and materials availability:”

Holds: README

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

Tools: NumPy (3 files), pandas (3 files), Matplotlib (2 files), PyBIDS (2 files), scikit-learn (2 files), SciPy (2 files), seaborn (2 files), statsmodels (2 files), XGBoost (2 files), neuroCombat (1 file)

Availability: 1 check, the latest on 7 November 2026: the link answers

7 November 2026: the link answers

4 files

SA_AD_classification.py, Python, 734 lines, 1 match
arc_harmonization.py, Python, 76 lines
clustering.py, Python, 1,691 lines, 2 matches
README.md, Text, 3 lines

At the source: github.com/zhangyubrain/Alzheimers-resilient-connectome

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;
3 scripts, each with its path and the digest of its content;
3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

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

Data

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

Data, code, and materials availability

All data needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. Requests for the data should be submitted to different public database portals. The OASIS-3 dataset is publicly available at www.oasis-brains.org/. The ADNI dataset is publicly available at https://adni.loni.usc.edu/. The HABS dataset is publicly available at https://habs.mgh.harvard.edu/. The code can be found at https://zenodo.org/records/20073656 and https://github.com/zhangyubrain/Alzheimers-resilient-connectome. This study did not generate new materials.

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

Versions

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

Version 1, 7 November 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 8 funders, 64 references.

Cite

This paper

Zhao, K., Xie, H., Fonzo, G. A., Carlisle, N. B., Jacobs, T., Osorio, R. S., Church, A., Lin, F. V., Zhang, Y., & Alzheimer’s Disease Neuroimaging Initiative. (2026). A functional connectome identifies vulnerability to cognitive decline and Alzheimer’s disease progression. Science advances, 12(45), eaed9351.

BibTeX

@article{zhao2026functional, author = {Zhao, Kanhao and Xie, Hua and Fonzo, Gregory A. and Carlisle, Nancy B. and Jacobs, Tovia and Osorio, Ricardo S. and Church, Arpana and Lin, Feng Vankee and Zhang, Yu and {Alzheimer’s Disease Neuroimaging Initiative}}, title = {{A functional connectome identifies vulnerability to cognitive decline and Alzheimer’s disease progression}}, journal = {Science advances}, year = {2026}, month = nov, volume = {12}, number = {45}, pages = {eaed9351}, publisher = {American Association for the Advancement of Science}, issn = {2375-2548}, pmcid = {PMC13733855} }

RIS

TY - JOUR AU - Zhao, Kanhao AU - Xie, Hua AU - Fonzo, Gregory A. AU - Carlisle, Nancy B. AU - Jacobs, Tovia AU - Osorio, Ricardo S. AU - Church, Arpana AU - Lin, Feng Vankee AU - Zhang, Yu AU - Alzheimer’s Disease Neuroimaging Initiative TI - A functional connectome identifies vulnerability to cognitive decline and Alzheimer’s disease progression T2 - Science advances J2 - Sci Adv PY - 2026 DA - 2026/11/06 VL - 12 IS - 45 SP - eaed9351 SN - 2375-2548 PB - American Association for the Advancement of Science LA - en ER -

CSL-JSON

{ "id": "pmcid:PMC13733855", "type": "article-journal", "title": "A functional connectome identifies vulnerability to cognitive decline and Alzheimer’s disease progression", "container-title": "Science advances", "author": [ { "family": "Zhao", "given": "Kanhao" }, { "family": "Xie", "given": "Hua" }, { "family": "Fonzo", "given": "Gregory A." }, { "family": "Carlisle", "given": "Nancy B." }, { "family": "Jacobs", "given": "Tovia" }, { "family": "Osorio", "given": "Ricardo S." }, { "family": "Church", "given": "Arpana" }, { "family": "Lin", "given": "Feng Vankee" }, { "family": "Zhang", "given": "Yu" }, { "literal": "Alzheimer’s Disease Neuroimaging Initiative" } ], "container-title-short": "Sci Adv", "volume": "12", "issue": "45", "page": "eaed9351", "PMCID": "PMC13733855", "ISSN": "2375-2548", "publisher": "American Association for the Advancement of Science", "language": "en", "issued": { "date-parts": [ [ 2026, 11, 6 ] ] } }

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.1001/jamanetworkopen.2024.20479 [code]: Dementia Subtypes Defined Through Neuropsychiatric Symptom-Associated Brain Connectivity Patterns.
Journal: JAMA network open
In common: PyBIDS, statsmodels, seaborn, 5 other tools, Alzheimer's / dementia, 4 references, 3 authors
[2] doi:10.1038/s41380-025-02974-6 [code]: Deep graph learning of multimodal brain networks defines treatment-predictive signatures in major depression.
Journal: Molecular psychiatry
In common: 5 references, 3 authors
[3] doi:10.1093/geront/gnaf277 [code]: What characterizes the exceptional cognition of superagers? A systematic review of multidomain biomarkers of successful cognitive aging.
Journal: The Gerontologist
In common: 9 references
[4] doi:10.1016/j.xcrm.2025.102151 [code]: Contrastive functional connectivity defines neurophysiology-informed symptom dimensions in major depression.
Journal: Cell reports. Medicine
In common: scikit-learn, pandas, SciPy, 2 other tools, 2 references, 2 authors
[5] doi:10.1038/s41398-022-02134-2: Transdiagnostic connectome signatures from resting-state fMRI predict individual-level intellectual capacity.
Journal: Translational psychiatry
In common: 3 references, 2 authors
[6] doi:10.1002/hbm.26453 [code]: Robust hierarchically organized whole-brain patterns of dysconnectivity in schizophrenia spectrum disorders observed after personalized intrinsic network topography.
Journal: Human brain mapping
In common: neuroCombat, PyBIDS, seaborn, 4 other tools, 5 references
[7] doi:10.1093/braincomms/fcae205 [code]: SuperAging functional connectomics from resting-state functional MRI.
Journal: Brain communications
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, Alzheimer's / dementia, 7 references
[8] doi:10.1016/j.dcn.2023.101337 [code]: A multi-sample evaluation of the measurement structure and function of the modified monetary incentive delay task in adolescents.
Journal: Developmental cognitive neuroscience
In common: PyBIDS, statsmodels, seaborn, 5 other tools, 5 references
[9] doi:10.1016/j.neurobiolaging.2022.06.008 [code]: Predicting future cognitive decline from non-brain and multimodal brain imaging data in healthy and pathological aging.
Journal: Neurobiology of aging
In common: PyBIDS, seaborn, scikit-learn, 4 other tools, Alzheimer's / dementia, 5 references
[10] doi:10.1007/s00234-026-04103-8 [code]: Enhanced detection of subtle cortical abnormalities in focal epilepsy using 7 T MRI surface-based models and graph neural networks.
Journal: Neuroradiology
In common: neuroCombat, PyBIDS, statsmodels, 6 other tools, 1 reference

Contribute

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

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

Add the badge to the code's README

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

Markdown

<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/repo/zenodo:20073656.svg)](https://openscicode.org/badge/repo/zenodo:20073656)
<!-- oscr-badge:end -->

reStructuredText

.. oscr-badge:start

.. image:: https://openscicode.org/badge/repo/zenodo:20073656.svg
   :target: https://openscicode.org/badge/repo/zenodo:20073656
   :alt: OSCR

.. oscr-badge:end

HTML

<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/repo/zenodo:20073656"><img src="https://openscicode.org/badge/repo/zenodo:20073656.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

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

zhangyubrain/Alzheimers-resilient-connectome: README.md on GitHub

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

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

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

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

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

Report a reproduction failure (signed in)

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

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