OSCR

Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,554 lines · 62 KB · no license

  1. # %% [markdown]
  2. # <a href="https://colab.research.google.com/github/dtabuena/TabuenaJangGrone/blob/main/_Clustering.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
  4. # %%
  5. !pip install openpyxl
  6. !pip install XlsxWriter
  7. !pip install pingouin
  8. !pip install CMH
  9. !pip install svgutils
  11. from IPython.display import clear_output
  12. import pandas as pd
  13. import os
  14. import numpy as np
  15. from sklearn.decomposition import PCA
  16. from matplotlib import pyplot as plt
  17. import matplotlib as mpl
  18. import seaborn as sns
  19. import matplotlib
  20. clear_output()
  22. from sklearn.cluster import KMeans
  23. import statsmodels.api as sm
  24. from statsmodels.formula.api import ols
  25. from statsmodels.stats.multicomp import pairwise_tukeyhsd
  26. import pingouin as pg
  27. import scipy
  28. # from pandas.compat.numpy import np_array_datetime64_compat
  29. from random import sample
  30. from tqdm import tqdm
  31. import os
  32. import shutil
  33. from google import colab
  35. from svgutils.compose import *
  37. clear_output()
  42. def restore_order(arg_ord):
  43. return np.argsort(arg_ord)
  45. # %%
  46. import urllib
  47. response = urllib.request.urlretrieve('https://raw.githubusercontent.com/dtabuena/Resources/main/Matplotlib_Config/Load_FS6.py','Load_FS6.py')
  48. %run 'Load_FS6.py'
  50. # %%
  52. def get_files(link):
  53. my_drop_folder = 'my_drop_folder'
  54. zipped_file_path = "/content/"+my_drop_folder + ".zip"
  55. unzipped_file_path = "/content/"+my_drop_folder
  56. if not( os.path.exists(zipped_file_path)):
  57. !wget -O $zipped_file_path $link # download with new name
  58. !echo A | unzip $zipped_file_path -d $unzipped_file_path
  59. file_list = [f for r,d,f in os.walk("/content/"+my_drop_folder )][0]
  60. return file_list
  62. def read_data_file(xl_to_analyze,dir='./Prelim_2023/Cleaned_DataSets/'):
  64. my_df_dict = pd.read_excel(dir+xl_to_analyze, engine='openpyxl',index_col=None,sheet_name=None)
  66. NEW_my_df_dict = {}
  67. for k in my_df_dict.keys():
  68. new_k = k
  69. if 'cre' in new_k:
  70. new_k = k.replace('Syn1-cre','/Syn1-Cre')
  71. NEW_my_df_dict[new_k] = my_df_dict[k]
  73. return NEW_my_df_dict, xl_to_analyze
  75. def dict_to_df(my_df_dict,AP_cut=60):
  76. '''convert dict of categories into a DF'''
  77. for k in my_df_dict.keys():
  78. df = my_df_dict[k]
  79. df['type'] = k
  80. df_list = [my_df_dict[k] for k in my_df_dict.keys()]
  81. full_df = pd.concat(df_list,ignore_index=True)
  83. '''drop low AP amps'''
  84. for r in full_df.index:
  85. if full_df.loc[r,'AP amp'] <AP_cut:
  86. full_df.at[r,'AP amp'] = np.nan
  87. return full_df
  89. def clean_data(labeling,full_df):
  91. numericals = [c for c in full_df.columns if c not in labeling ]
  93. full_df_norm = full_df.copy()
  94. mean_dict = {}
  95. sd_dict = {}
  97. if 'extras' not in full_df.columns:
  98. exclude_extras = False
  99. plot_extras = False
  101. for c in full_df_norm.columns:
  102. if c not in labeling:
  103. col_mean = np.mean(full_df_norm[c])
  104. col_sd = np.std(full_df_norm[c])
  105. mean_dict[c] = col_mean
  106. sd_dict[c] = col_sd
  107. full_df_norm[c] = (full_df_norm[c] - col_mean)/col_sd
  110. null_check = pd.isna(full_df_norm)
  111. incomplete =[]
  112. for r in full_df_norm.index:
  113. vals = np.array(null_check.loc[r,numericals])
  114. bad = any(vals)
  115. if bad:
  116. incomplete.append(r)
  119. complete = [c for c in full_df_norm.index if c not in incomplete]
  120. full_df_no_null = full_df_norm.loc[complete].copy()
  123. if exclude_extras:
  124. is_extra = list(full_df_no_null['extra']==1)
  125. extra_ind = full_df_no_null.index[is_extra]
  126. full_df_no_null.drop( index =extra_ind, inplace=True)
  129. data_X = full_df_no_null[numericals].to_numpy()
  130. data_for_fit = data_X
  131. numericals = [n.replace('(1)', '') for n in numericals]
  132. return numericals, data_for_fit, full_df_no_null
  134. def get_types(full_df_no_null):
  135. type_list = list(full_df_no_null['type'])
  136. uniq_types = list(set(type_list))
  137. young_types = [t for t in uniq_types if '-9' in t and 'Cre' not in t]
  138. young_types.sort()
  139. cre_types = [t for t in uniq_types if 'Cre' in t]
  140. cre_types.sort()
  141. old_types = [t for t in uniq_types if '-19' in t and 'Cre' not in t]
  142. old_types.sort()
  143. uniq_types = young_types + cre_types + old_types
  144. # print(uniq_types)
  145. num_t = len(uniq_types)
  146. type_num = [ uniq_types.index(t) for t in type_list]
  147. return type_num, num_t, uniq_types,type_list
  149. def cluster_full(data_for_fit,xl_to_analyze,labels):
  150. random_state = 42
  151. if 'CA1' in xl_to_analyze:
  152. random_state = 47
  153. if 'Type I ' in xl_to_analyze:
  154. random_state = 43
  155. if 'Type II ' in xl_to_analyze:
  156. random_state = 42
  159. rheo_ind = [i for i in range(len(labels)) if 'Rheo' in labels[i]][0]
  161. # kmeans = KMeans(n_clusters=2, random_state=random_state).fit(data_for_fit)
  162. needs_flip = True
  163. while needs_flip:
  164. random_state+=1
  165. print('random_state',random_state)
  166. kmeans = KMeans(n_clusters=2, random_state=random_state,n_init=10).fit(data_for_fit)
  167. centers = kmeans.cluster_centers_.T
  168. needs_flip = np.diff(centers[rheo_ind])<0
  169. print(np.diff(centers[rheo_ind]))
  171. # centers = kmeans.cluster_centers_.T
  172. distances = kmeans.transform(data_for_fit)
  173. rel_centers = (centers.T - np.mean(centers.T,axis = 0)).T
  176. labels = kmeans.labels_
  177. rel_counts = [np.mean( l == labels ) for l in range(2)]
  178. # new_cent_order = np.argsort(centers[:,0]).tolist()
  179. new_cent_order = np.arange(len(centers)).tolist()
  181. total_D = np.tile(np.sum(distances,axis=1),[2,1]).T
  182. closeness = np.log(total_D/distances)
  183. closeness = (closeness - np.mean(closeness,axis=0)) / (np.std(closeness,axis=0))
  184. distances_norm = (distances - np.mean(distances,axis=0)) / (np.std(distances,axis=0))
  185. score = closeness[:,1] - closeness[:,0]
  187. packaged = {'features':numericals,
  188. 'kmeans':kmeans,
  189. 'centers':centers,
  190. 'distances':distances,
  191. 'rel_centers':rel_centers,
  192. 'rel_counts':rel_counts,
  193. 'total_D':total_D,
  194. 'closeness':closeness,
  195. 'distances_norm':distances_norm,
  196. 'score':-score,
  197. 'data_for_fit':data_for_fit,
  198. 'new_cent_order':new_cent_order
  199. }
  201. return packaged
  205. def do_stats(type_list,score):
  206. is_E3 = ['E3' in t for t in type_list ]
  207. is_9mo = ['-19' not in t for t in type_list ]
  208. age_list = ['7-9' if b else '17-19' for b in is_9mo]
  209. genotype = ['apoE3-KI' if b else 'apoE4-KI' for b in is_E3]
  210. is_crePos = ['cre+' in t for t in type_list ]
  211. is_creNeg = ['cre-' in t for t in type_list ]
  212. not_cree = [not any( [is_crePos[i], is_creNeg[i]] ) for i in range(len(is_crePos))]
  213. e34_genotype = [genotype[i] for i in range(len(not_cree)) if not_cree[i]]
  214. e34_score = [score[i] for i in range(len(not_cree)) if not_cree[i]]
  215. e34_age = [str(age_list[i])+'mo' for i in range(len(not_cree)) if not_cree[i]]
  216. genotype_list_str = [t.replace(' 7-9mo', '').replace(' 17-19mo', '') for t in type_list]
  217. age_list_str = [str(a)+'mo' for a in age_list]
  219. df_for_34_2W_anova = pd.DataFrame({'genotype': genotype_list_str, 'Age': age_list_str,'T1_Score': score})
  222. ANOVA2_results = pg.anova(dv='T1_Score', between=['genotype', 'Age'], data=df_for_34_2W_anova, detailed=True)
  223. # display(ANOVA2_results)
  225. ANOVA2_results.to_csv(f'{folder}/2W_results Age-Geno.csv')
  226. df_for_34_2W_anova['combination'] = df_for_34_2W_anova.genotype + " " + df_for_34_2W_anova.Age
  227. m_comp = pairwise_tukeyhsd(endog=df_for_34_2W_anova['T1_Score'], groups=df_for_34_2W_anova['combination'], alpha=0.05)
  228. tukey_data = pd.DataFrame(data=m_comp._results_table.data[1:], columns = m_comp._results_table.data[0])
  229. tukey_data = tukey_data.sort_values('p-adj',ascending=True)
  230. # display(tukey_data)
  232. tukey_data.to_csv(f'{folder}/pairwise Age-Genotype.csv')
  233. df_for_34_2W_anova.to_csv(f'{folder}/Scores.csv')
  234. return df_for_34_2W_anova, tukey_data
  237. def save_dl(folder):
  238. root_list = []
  239. f_list = []
  241. fold_name = folder
  242. fold_name = fold_name.replace(' ','_')
  243. zip_name= fold_name+'.zip'
  245. print(folder)
  246. !zip -r $zip_name $folder
  247. colab.files.download(zip_name)
  248. return None
  252. # %%
  253. def plot_summary(num_t,numericals,type_list,uniq_types,packaged,df_for_34_2W_anova,tukey_data,folder, clust_names = ['Hyper', 'Normal'],fig_opt = {'dpi': 300, 'format':'png','bbox_inches':None}):
  255. cmap = plt.cm.cool(np.linspace(0,1,2))*np.array([0.6]*3+[1])
  256. dark_cool = mpl.colors.ListedColormap(cmap)
  257. dpi=300
  259. kmeans = packaged['kmeans']
  260. centers = packaged['centers']
  261. distances = packaged['distances']
  262. rel_centers = packaged['rel_centers']
  263. rel_counts = packaged['rel_counts']
  264. total_D = packaged['total_D']
  265. closeness = packaged['closeness']
  266. distances_norm = packaged['distances_norm']
  267. score = packaged['score']
  268. data_for_fit = packaged['data_for_fit']
  269. cell_type = packaged['cell_type']
  270. type_num = packaged['type_num']
  272. colors = plt.cm.viridis(np.linspace(0,1,num_t))
  273. alpha_1 = 1
  274. alpha_2 = .8
  275. colors = [[4/255, 51/255, 255/255, alpha_1], [4/255, 51/255, 255/255, alpha_2],
  276. [255/255, 38/255, 0/255, alpha_1], [255/255, 38/255, 0/255, alpha_2],
  277. [103/255, 35/255, 2/255, alpha_1], [103/255, 35/255, 2/255, alpha_2],
  278. ]
  281. list_colors = [colors[i] for i in type_num ]
  282. color_dc={'apoE3-KI 7-9mo': colors[0],
  283. 'apoE3-KI 17-19mo': colors[1],
  284. 'apoE4-KI 7-9mo':colors[2],
  285. 'apoE4-KI 17-19mo': colors[3],
  286. 'fE4/Syn1-Cre+ 7-9mo': colors[4],
  287. 'fE4/Syn1-Cre- 7-9mo': colors[5],}
  290. ################# raw_data_map '''
  291. counts = [np.sum( [tt==t for tt in type_list]) for t in uniq_types ]
  292. counts[-1] = counts[-1] + np.sum(counts)*0.075
  293. fig_norm_vals, ax = plt.subplots(1,num_t+1,figsize = [7.25,1.6],gridspec_kw={'width_ratios': counts+[.75]},layout="constrained",dpi=dpi)
  294. data_X_min = np.min(data_for_fit)
  295. data_X_max = np.max(data_for_fit)
  296. extrema = np.array([-1,1]) * np.max(np.abs([data_X_min,data_X_max]))
  297. data_X_min = extrema[0]
  298. data_X_max = extrema[1]
  300. for ti in range(num_t):
  301. t = uniq_types[ti]
  302. t_colr = colors[ti]
  303. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  304. c = ax[ti].pcolor(data_for_fit.T[:,is_type], vmin=data_X_min*1.2, vmax=data_X_max*1.2)
  305. ax[ti].set_yticks(np.arange(len(numericals))+0.5)
  306. if ti ==0 : ax[ti].set_yticklabels(numericals,rotation=0)
  307. else: ax[ti].get_yaxis().set_visible(False)
  308. if ti==0: ax[ti].set_xlabel('Cell (#)')
  309. ax[ti].title.set_text(t.replace(' ','\n'))
  310. start = min(is_type)+3 / len(is_type)
  311. stop = max(is_type)-3 / len(is_type)
  312. fig_norm_vals.suptitle(cell_type+ ' Normalized Parameters',y=1.05,)# fontsize=6
  314. plt.colorbar(c, cax=ax[ti+1],label='z-scored'+'\n'+'values')
  315. plt.show()
  316. # plt.tight_layout()
  319. # plt.rcParams.update({'font.size': FS,'font.family': 'arial'}) #14
  322. ################# Centers
  323. new_cent_order = np.argsort(centers[:,1]).tolist()[::-1]
  324. fig_centers, ax = plt.subplots(1,1,figsize = [3,2],dpi=dpi,layout="constrained")
  325. ax.set_prop_cycle(color=[cmap[0], cmap[1]])
  326. ax.plot(centers[new_cent_order,:],'-o',linewidth=.5,markersize=2)
  327. ax.legend(clust_names,loc='lower center',frameon=False)
  328. sorted_numericals = [numericals[i] for i in new_cent_order]
  329. ax.set_xticks(np.arange(len(numericals)))
  330. ax.set_xticklabels(sorted_numericals,rotation=45,ha='right', rotation_mode='anchor')
  331. ax.set_ylabel('z-scored'+' '+'values')
  332. ax.axhline(0,color='k',linewidth=1)
  333. ax.title.set_text(cell_type +' '+'Cluster Centers')
  334. # plt.tight_layout()
  336. ################# Split Distance heat map
  337. counts = [np.sum( [tt==t for tt in type_list]) for t in uniq_types ]
  338. print(counts)
  339. b_dim = 2
  340. fig_distances_split, ax = plt.subplots(1,num_t+1,figsize = [7.25,.9],gridspec_kw={'width_ratios': counts+[.75]},layout="constrained",dpi=dpi) ## NEw
  341. dist_min = np.min(distances_norm)
  342. dist_max = np.max(distances_norm)
  344. extrema = np.array([-1,1]) * np.max(np.abs([dist_min,dist_max]))
  345. dist_min = extrema[0]
  346. dist_max = extrema[1]
  348. # plt.rcParams.update({'font.size': FS,'font.family': 'arial'}) #14
  349. for ti in range(num_t):
  350. t = uniq_types[ti]
  351. t_colr = colors[ti]
  352. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  353. c = ax[ti].pcolor(distances_norm.T[:,is_type], vmin=dist_min*1.2, vmax=dist_max*1.2)
  354. ax[ti].set_yticks(np.arange(len(clust_names))+0.5)
  355. clust_names_nl = [c.replace(' ',' ') for c in clust_names]
  356. if ti ==0 : ax[ti].set_yticklabels(clust_names_nl,rotation=0)
  357. else: ax[ti].get_yaxis().set_visible(False)
  358. if ti==0: ax[ti].set_xlabel('Cell (#)')
  359. # fig_distances_split.tight_layout()
  360. fig_distances_split.colorbar(c, cax=ax[ti+1],label='z-scored'+'\n'+'distance')
  361. for ti in range(num_t):
  362. t = uniq_types[ti]
  363. ax[ti].text(0.5, 1, t.replace(' ','\n'),transform=ax[ti].transAxes,ha='center',rotation=0,va='bottom') # ,fontsize=14
  364. # fig_distances_split.tight_layout()
  365. plt.show()
  368. ############# Combined Scatter
  369. # plt.rcParams.update({'font.size': FS,'font.family': 'arial'}) #14
  370. fig_comb_scatter, ax = plt.subplots(1,figsize = [7.25/5,7.25/5],dpi=dpi,layout="constrained")
  371. for ti in range(num_t):
  372. t = uniq_types[ti]
  373. t_colr = color_dc[t]
  374. print(t,t_colr)
  375. t_colr_edge = t_colr[:-1] + [alpha_1]
  376. t_colr_edge = None
  377. t_colr = 'k'
  378. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  379. ax.scatter(distances_norm[is_type,0],distances_norm[is_type,1],color=t_colr, label=t,linewidth=1,edgecolor=t_colr_edge,s=2,zorder=1)
  380. # ax.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
  381. ax.set_xlabel('Dist. from '+clust_names[0]+' Center \n (z-scored)')
  382. ax.set_ylabel('Dist. from '+clust_names[1]+' Center \n (z-scored)')
  383. ax.plot([0, 1], [0, 1], transform=ax.transAxes,color='k',linewidth=1)
  384. ax.text(0.98, 0.02, ''+clust_names[1],transform=ax.transAxes,ha='right',va='bottom') # ,fontsize=12
  385. ax.text(0.02, 0.98, ''+clust_names[0],transform=ax.transAxes,ha='left',va='top') # ,fontsize=12
  386. comb_xlim = ax.get_xlim()
  387. comb_ylim = ax.get_ylim()
  388. common_lim = np.array([-1,1]) * np.max(np.abs(list(comb_xlim)+list(comb_ylim)))
  389. common_lim= [np.min([comb_xlim[0],comb_ylim[0]]), np.max([comb_xlim[1],comb_ylim[1]]) ]
  390. ax.set_xlim(common_lim)
  391. ax.set_ylim(common_lim)
  393. # residuals = np.diff(distances_norm,axis=1)
  394. x=distances_norm[:,0]
  395. y=distances_norm[:,1]
  396. for i in range(len(x)):
  397. if x[i]>y[i]:
  398. c=[.6, 0, .6]
  399. else :
  400. c=[0, .6, .6]
  401. plt.plot([x[i], x[i]], [y[i], x[i]], color=c, linewidth=0.5,zorder=-1)
  403. #### separate Scatter #####
  404. b_dim = 3
  406. print(common_lim)
  407. common_lim =common_lim+np.array([0,2])
  409. """ for compact scatter Layout """
  410. if 'CA3' in cell_type: fig_sep_scatter, axs = plt.subplots(3,int(num_t/3),figsize = [7.25/4,7.25/4*1.5*1.2],dpi=dpi,layout="constrained") #
  411. else: fig_sep_scatter, axs = plt.subplots(2,int(num_t/2),figsize =[7.25/4,7.25/4*1.2],dpi=dpi,layout="constrained") # [int(num_t/2)*b_dim,b_dim*2]
  413. axs = axs.T.flatten()
  414. if 'CA3' in cell_type:
  415. axs = axs[np.arange(6)[[0,1,2,5,3,4]]]
  417. for ti in np.arange(num_t):
  418. t = uniq_types[ti]
  419. t_colr = color_dc[t]
  420. t_colr_edge = t_colr[:-1] + [alpha_1]
  421. t_colr_edge = None
  422. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  423. axs[ti].scatter(distances_norm[is_type,0],distances_norm[is_type,1],color=t_colr, label=t,linewidth=1,edgecolor=t_colr_edge,s=2)
  424. axs[ti].plot([0, 1], [0, 1], transform=axs[ti].transAxes,color=[0,0,0,1],linewidth=1)
  425. axs[ti].title.set_text(t.replace(' ','\n').replace('apo',""))
  426. axs[ti].set_xlim(common_lim)
  427. axs[ti].set_ylim(common_lim)
  429. axs[ti].text(0.98, 0.02, ''+clust_names[1],transform=axs[ti].transAxes,ha='right',va='bottom')
  430. axs[ti].text(0.02, 0.98, ''+clust_names[0],transform=axs[ti].transAxes,ha='left',va='top')
  433. fig_sep_scatter.subplots_adjust(wspace=0.3, hspace=0.3)
  434. fig_sep_scatter.supxlabel(f'{clust_names[0]} Distance', x=0.60,y=.05)
  435. fig_sep_scatter.supylabel(f'{clust_names[1]} Distance', x=0.050)
  436. plt.tight_layout()
  441. ######## ''' Violins of Residual Plots '''
  442. if 'CA3' in cell_type: fig_violin_residuals, ax = plt.subplots(1,figsize=np.array([7.25*.3,2.25]),dpi=dpi,layout="constrained")
  443. else: fig_violin_residuals, ax = plt.subplots(1,figsize=np.array([7.25*.7/3,2.25]),dpi=dpi,layout="constrained")
  445. comb_list = list(set(df_for_34_2W_anova['combination']))
  447. young_types = [t for t in comb_list if '-9' in t and 'Cre' not in t]
  448. young_types.sort()
  449. cre_types = [t for t in comb_list if 'Cre' in t]
  450. cre_types.sort()
  451. old_types = [t for t in comb_list if '-19' in t and 'Cre' not in t]
  452. old_types.sort()
  453. comb_list = young_types + cre_types + old_types
  455. sns_dict={}
  459. for ti in range(len(comb_list)):
  460. type_i = comb_list[ti]
  461. match_type = df_for_34_2W_anova['combination']==type_i
  462. y_vals = list(np.array(df_for_34_2W_anova[match_type]['T1_Score']))
  463. x_vals = np.ones_like(y_vals)*ti + int(np.floor(ti/2))
  464. pos = list(np.ones_like(y_vals)*ti)
  465. t_colr = colors[ti]
  466. t_colr = color_dc[type_i]
  467. sns_dict[type_i] =y_vals
  469. sns_df = pd.DataFrame.from_dict(sns_dict, orient='index').T
  471. colors_sns = [color_dc[t] for t in sns_df.columns]
  472. colors_sns = [ list(np.array(c[:-1]) + (1 - np.array(c[:-1]))* (1-c[-1])) + list(c[-1:]) for c in colors_sns]
  474. sns.swarmplot(ax=ax,data=sns_df,palette=colors_sns, label=t,linewidth=.05,edgecolor='k',size=2) # color='k'
  475. sns.violinplot(ax=ax,data=sns_df,palette=colors_sns,scale ='area', label=t, linewidth=1,inner=None)
  476. # Manually adjust the violins
  477. violins = [v for v in ax.collections if isinstance(v, matplotlib.collections.PolyCollection)] # Get all violins
  478. for i, violin in enumerate(violins):
  479. violin.set_facecolor([0.9, 0.9, 0.9]) # Set the fill color to light gray
  480. violin.set_edgecolor(colors_sns[i % len(colors_sns)]) # Set edge color using colors_sns
  481. violin.set_linewidth(1) # Set the linewidth to make the edges more visible
  483. ax.text(0, 1, clust_names[0],transform=ax.transAxes,ha='right',va='bottom',color=[0, .6, .6])
  484. ax.text(0, 0, clust_names[1],transform=ax.transAxes,ha='right',va='top',color=[.6, 0, .6])
  486. ax.set_ylabel('Residual Value')
  487. ax.set_ylim(-7.5,12.5)
  488. ax.axhline(0,color='k',linewidth=.5)
  489. ax_labels = sns_df.columns
  490. ax.set_xticks(ticks=range(len(ax_labels)))
  491. ax_labels_r = [c.replace(' ','\n') for c in ax_labels]
  492. ax_labels_r = [c.replace('\nsyn',' syn') for c in ax_labels_r]
  493. ax_labels_r = [c.split('\n')[0] for c in ax_labels_r]
  494. ax_labels_r = [c.replace('fE4/Syn1-Cre', 'Cre') for c in ax_labels_r]
  495. ax_labels_r = [c.replace('apo', '') for c in ax_labels_r]
  496. # ax_labels_r = ['\n'+ax_labels_r[i] if i%2 == 0 else ax_labels_r[i] for i in range(len(ax_labels_r)) ]
  499. # ax.set_xticklabels(ax_labels_r,rotation=60,ha='right', va='top', rotation_mode='anchor')
  500. ax.set_xticklabels(ax_labels_r,rotation=0,ha='center', va='top', rotation_mode='anchor')
  501. ax.title.set_text(cell_type) #+' Cluster Membership'
  503. div_x = np.sum( ['7-9' in l for l in ax_labels])-0.5
  504. segment_x = [div_x, div_x]
  505. segment_y = ax.get_ylim()*np.array([1,1])
  506. ll = plt.plot(segment_x,segment_y,':k',linewidth=.5) # returns a sequence of line objects
  507. ll[0].set_clip_on(False)
  508. month_offset = 0 # 10
  509. plt.text(div_x/2-0.25, segment_y[0]-month_offset, '7-9mo',ha='center',va='bottom')
  510. plt.text(div_x+1, segment_y[0]-month_offset, '17-19mo', ha='center',va='bottom')
  511. print('sns_dict',sns_dict)
  513. cre_ind = [i for i,k in enumerate(sns_dict.keys()) if 'Cre' in k]
  514. if len(cre_ind)>0:
  515. x = np.mean(cre_ind)
  516. plt.text(x, segment_y[0]-4, 'fE4/Syn-Cre',ha='center',va='bottom')
  517. else:
  518. plt.text(1, segment_y[0]-4, 'PLACE HOLDER' ,ha='center',va='bottom', color = 'white')
  520. base_offset = 1.5
  521. y0 = np.max(df_for_34_2W_anova['T1_Score'])*base_offset
  522. print('y0',y0)
  523. labels = comb_list
  524. true_tuckey = tukey_data[tukey_data['reject']==True]
  526. astk_range = [0.05,0.01,0.001,0.0001]
  527. for r in true_tuckey.index:
  528. comps = ['apoE3-KI 7-9mo']
  529. if true_tuckey.loc[r,'group1'] in comps or true_tuckey.loc[r,'group2'] in comps:
  530. start_l = true_tuckey.loc[r,'group1']
  531. stop_l = true_tuckey.loc[r,'group2']
  532. start_x = labels.index(start_l)
  533. stop_x = labels.index(stop_l)
  534. y0 = np.max(df_for_34_2W_anova['T1_Score'])*base_offset - 1
  535. y_offset = abs((stop_x-start_x))
  536. y = y0+ y_offset * 1.25 + np.min([stop_x,start_x])*.25
  538. half_dx = abs(stop_x-start_x)/2
  539. ax.errorbar( np.mean([stop_x,start_x]), y,yerr=0, xerr=half_dx,color='k',capsize=.75,linewidth=0.75,capthick=0.75)
  541. cent =np.mean([start_x,stop_x])
  542. num_ask = np.sum([true_tuckey.loc[r,'p-adj'] < t for t in astk_range] )
  544. plt.text(cent,y, '*' * num_ask,ha='center',va='center',fontsize = 2+plt.rcParams['font.size'])
  547. if 'Syn' in true_tuckey.loc[r,'group1'] and 'Syn' in true_tuckey.loc[r,'group2']:
  548. start_l = true_tuckey.loc[r,'group1']
  549. stop_l = true_tuckey.loc[r,'group2']
  550. start_x = labels.index(start_l)
  551. stop_x = labels.index(stop_l)
  552. y0 = np.max(df_for_34_2W_anova['T1_Score'])*base_offset
  553. y_offset = abs((stop_x-start_x))
  554. y = y0+ y_offset * .75
  556. half_dx = abs(stop_x-start_x)/2
  557. ax.errorbar( np.mean([stop_x,start_x]), y,yerr=0, xerr=half_dx,color='k',capsize=.75,linewidth=0.75,capthick=0.75)
  559. cent =np.mean([start_x,stop_x])
  560. num_ask = np.sum([true_tuckey.loc[r,'p-adj'] < t for t in astk_range] )
  561. print(true_tuckey.loc[r,'p-adj'],num_ask)
  562. plt.text(cent,y, '*' * num_ask,ha='center',va='center',fontsize = 2+plt.rcParams['font.size'])
  564. plt.tight_layout()
  567. fig_norm_vals.savefig(f'{folder}/Normalized_Vals_split.'+fig_opt['format'],**fig_opt)
  568. fig_violin_residuals.savefig(f'{folder}/fig_violin_residuals.'+fig_opt['format'],**fig_opt)
  569. fig_distances_split.savefig(f'{folder}/fig_distances_split.'+fig_opt['format'],**fig_opt)
  570. fig_sep_scatter.savefig(f'{folder}/GenoType_dist_Separate.'+fig_opt['format'],**fig_opt)
  571. fig_comb_scatter.savefig(f'{folder}/GenoType_dist.'+fig_opt['format'],**fig_opt)
  573. return None
  580. # %%
  581. import matplotlib.pyplot as plt
  582. import numpy as np
  583. from matplotlib.patches import FancyArrow
  585. def fig_cartoon_x(fig_opt={'dpi': 300, 'format': 'svg', 'bbox_inches': None}):
  586. FS = 8
  587. A = np.random.randn(16, 80)
  588. B = np.random.randn(2, 80)
  589. B_t = np.random.randn(16, 2)
  591. fig_cartoon, ax = plt.subplots(3, 4, figsize=[7.25 * .85, .5], width_ratios=[6, 3, .2, 6], height_ratios=[4, 2, 4], dpi=300, layout="constrained")
  593. # Hide the original axes
  594. for a in ax.reshape(-1):
  595. a.set_visible(False)
  597. # Get the gridspec to modify the layout
  598. gs = ax[0, 0].get_gridspec()
  600. # Create the large axis for the heatmap A
  601. big_ax_1 = fig_cartoon.add_subplot(gs[0:, 0])
  602. ch = big_ax_1.pcolorfast(A)
  603. cbar = plt.colorbar(ch)
  604. cbar.set_ticks([])
  605. cbar.set_label('z-score')
  606. big_ax_1.set_xlabel('Pooled CA3 Cells', fontsize=FS)
  607. big_ax_1.set_ylabel('16 Parameters', rotation=0, ha='right', va='center', fontsize=FS)
  608. big_ax_1.xaxis.set_label_position('top')
  609. big_ax_1.set_xticks([])
  610. big_ax_1.set_yticks([])
  612. # Create an axis for B (second row, fourth column)
  613. ax[1, 3].set_visible(True)
  614. ax[1, 3].pcolorfast(B)
  615. ax[1, 3].set_ylabel('k', rotation=0, ha='right', va='center', fontsize=FS)
  616. ax[1, 3].set_title('Cell Distances', rotation=0, ha='center', va='center', fontsize=FS)
  617. ax[1, 3].set_xticks([])
  618. ax[1, 3].set_yticks([])
  620. # Create the third large axis for B_t
  621. big_ax_2 = fig_cartoon.add_subplot(gs[0:, 2])
  622. big_ax_2.pcolorfast(B_t)
  623. big_ax_2.set_xlabel('k', rotation=0, ha='center', va='center', fontsize=FS)
  624. big_ax_2.set_ylabel('Centroid\n Parameters (16)', rotation=90, ha='center', va='bottom', fontsize=FS)
  625. big_ax_2.set_xticks([])
  626. big_ax_2.set_yticks([])
  628. # Create a new axis for the blocky arrow, spanning only the second column
  629. arrow_ax = fig_cartoon.add_subplot(gs[:, 1]) # this spans the whole second column
  630. arrow_ax.set_xticks([])
  631. arrow_ax.set_yticks([])
  632. arrow_ax.set_frame_on(False) # Hide the axis frame
  633. arrow_ax.set_xlim(-.0,1.2)
  635. # Add the blocky arrow on the new axis
  636. x_start, y_start = 0, 0.5 # Arrow starting position in axes coordinates
  637. dx, dy = 0.6, 0 # Arrow length (horizontal)
  638. arrow = FancyArrow(x_start, y_start, dx, dy, width=.35, head_width=0.75, head_length=0.15, color='k', transform=arrow_ax.transAxes)
  639. arrow_ax.add_patch(arrow)
  641. # Add label for the arrow
  642. arrow_ax.annotate('k-means', xy=(x_start + dx / 2, y_start), ha='center', va='center', xycoords='axes fraction', color='w', fontsize=FS)
  644. # Save the figure
  645. fig_cartoon.savefig('./fig_cartoon_x.' + fig_opt['format'], **fig_opt)
  647. return None
  652. # %%
  653. def cross_val(data_for_fit,numericals,new_cent_order,folder,k=2,
  654. num_iter = 1000,SubSampleFrac = .8,cell_type='',fig_opt = {'dpi': 300, 'format':'svg','bbox_inches':None},ds_iter=10):
  656. cval_res = {}
  657. center_list = []
  658. label_list = []
  659. sub_sample_ind_list = []
  660. full_label_list = []
  661. pred_all = True
  662. # sorted_numericals = [numericals[i] for i in new_cent_order]
  664. FS=6
  665. dpi=300
  667. # set_font_all(FS)
  669. cmap = plt.cm.cool(np.linspace(0,1,k))*np.array([0.6]*3+[1])
  670. dark_cool = mpl.colors.ListedColormap(cmap)
  672. for ni in tqdm( range(num_iter) ):
  673. sub_sample_ind = []
  674. for t in uniq_types:
  675. t_inds = [i for i in range(len(type_list)) if t in type_list[i]]
  676. # print(t_inds)
  677. sub_sample_ind.extend( sample(t_inds, int(len(t_inds)*SubSampleFrac) ) )
  678. kmeans_iter = KMeans(n_clusters=k, random_state=42,n_init=10).fit(data_for_fit[sub_sample_ind][:,new_cent_order])
  679. centers_iter = kmeans_iter.cluster_centers_.T
  680. centers_iter = centers_iter - np.mean(centers_iter,axis=0)
  681. labels_iter = kmeans_iter.labels_
  682. if pred_all:
  683. labels_iter = kmeans_iter.predict(data_for_fit)
  685. if ni>0:
  686. corco = np.corrcoef(centers_iter.T,center_list[0].T)
  687. corco = corco[0:k,k:]
  688. ind = np.argmax(corco,axis=0)
  689. centers_iter = centers_iter[:,ind]
  690. labels_iter = [ind[l] for l in labels_iter]
  692. if not pred_all:
  693. full_label_iter = np.empty((data_for_fit.shape[0],))
  694. full_label_iter[:] = np.nan
  695. full_label_iter[sub_sample_ind] = labels_iter
  696. else:
  697. full_label_iter = labels_iter
  699. center_list.append(centers_iter)
  700. label_list.append(labels_iter)
  701. sub_sample_ind_list.append(sub_sample_ind_list)
  702. full_label_list.append(full_label_iter)
  704. center_stack = np.stack(center_list)
  705. mean_center = np.mean(center_stack,0)
  706. nintey_p_center = np.percentile(center_stack, [.5, 99.5], axis=0)
  708. cval_res['mean_center'] = mean_center
  709. cval_res['confidence'] = nintey_p_center
  712. plot_order = np.argsort(np.diff(mean_center,1).T)[0]
  713. numericals=np.array(numericals)
  715. ####### PLOT CENTERS
  716. plt.rcParams.update({'font.size': FS,'font.family': 'arial'}) #14
  717. fig_cv_cent, ax = plt.subplots(1,figsize = (7.25/4,1.8),dpi=dpi,layout="constrained")
  719. handle_list = []
  720. if k == 2:
  721. clust_names = [ 'Hyper', 'Normal']
  722. rheo_index_in_center = [i for i in range(len(numericals)) if 'Rheo' in numericals[i]]
  723. rheo_dif = np.diff(mean_center[rheo_index_in_center,:])[0]
  724. if rheo_dif<0:
  725. clust_names=clust_names[::-1]
  726. cmap=np.flipud(cmap)
  728. for ki in range(k):
  729. if k==2: ax.plot(mean_center[plot_order,ki].T,'-o',color=cmap[ki],linewidth=1,label = clust_names[ki],markersize=2)
  730. else: ax.plot(mean_center[plot_order,ki].T,'-o',color=cmap[ki],linewidth=.5,markersize=2)
  731. ax.plot(nintey_p_center[0,plot_order,ki].T,':',color=cmap[ki],linewidth=.5)
  732. ax.plot(nintey_p_center[1,plot_order,ki].T,':',color=cmap[ki],linewidth=.5)
  734. ax.set_xticks(np.arange(len(numericals)))
  735. ax.set_xticklabels(numericals[plot_order],rotation=45,ha='right',va='center', rotation_mode='anchor',fontsize=plt.rcParams['font.size']) # ,fontsize=14
  736. ax.axhline(0, color='k',linewidth=.5)
  738. if k == 2: ax.legend(loc='upper center',frameon=False,fontsize=plt.rcParams['font.size'])
  739. ax.set_ylabel('z-scored'+' '+'values')
  742. full_label_stack = np.stack(full_label_list)
  744. pseudo_truth = scipy.stats.mode(full_label_stack,axis=0, nan_policy='omit')[0]
  745. is_match = (full_label_stack == pseudo_truth)*1. # convert to float so can accept nans
  746. is_match[np.isnan(full_label_stack)]=np.nan
  747. consistency_iter = np.nanmean(is_match, axis = 1)
  748. mean_consistency_iter = np.nanmean(consistency_iter)
  749. std_consistency_iter = np.nanstd(consistency_iter)
  750. consistency_cell = np.nanmean(is_match, axis = 0)
  751. mean_consistency_cell = np.nanmean(consistency_cell)
  752. sd_consistency_cell = np.nanstd(consistency_cell)
  753. consistency_report = 'label consistency: ' + str(round(mean_consistency_cell*100,1)) + '% +- ' + str(round(sd_consistency_cell*100,2)) + ' (chance: ' + str(round(1/k*100,1))+')'
  754. ax.title.set_text(cell_type)
  757. cval_res['mean_consistency_iter'] = mean_consistency_iter
  758. cval_res['std_consistency_iter'] = std_consistency_iter
  759. cval_res['consistency_cell'] = consistency_cell
  760. cval_res['mean_consistency_cell'] = mean_consistency_cell
  762. plt.tight_layout()
  763. fig_cv_cent.tight_layout()
  764. fig_cv_cent.patch.set_facecolor('white')
  766. plt.show()
  767. fig_cv_cent.savefig(f'{folder}/CrossVal_Data k={k} Centers.'+fig_opt['format'],**fig_opt)
  770. ################## Plot iteration results
  771. scale_factor= 1
  772. plt.rcParams.update({'font.size': FS*scale_factor, 'font.family': 'arial'})
  773. fig_cv_prob, ax = plt.subplots(2,2,figsize=np.array([7.25/4,1.8]),
  774. gridspec_kw={'height_ratios': [3, 1],'width_ratios': [4, 1.1]} ,dpi=dpi,layout="constrained")
  776. ax = ax.T.flatten()
  778. iters_ds = np.arange(0, num_iter + 1, ds_iter)
  779. cell_ids = np.arange(0, full_label_stack.shape[1] + 1) # **Use shape[1] for cells**
  780. pcol = ax[0].pcolorfast(cell_ids, iters_ds, full_label_stack[::ds_iter,:], cmap=dark_cool) # **Swap arguments**
  783. ax[0].set_xlabel('Cell (#)\n',va='bottom')
  784. ax[0].set_ylabel('Iteration (#)\n',va='bottom',rotation = 90)
  785. ax[0].xaxis.set_label_position('top')
  786. # ax[0].set_xticks([])
  787. ax[0].xaxis.tick_top()
  789. ax[1].scatter(np.arange(len(consistency_cell)),consistency_cell*100,color='k',s=1,edgecolor='k',marker='.')
  790. ax[1].axhline(1/k*100,color='k',linestyle=':',linewidth=.25)
  791. ax[1].set_ylabel('Cell\nConsistency (%)',va='bottom',rotation = 90)
  792. ax[1].set_ylim((0,102))
  793. ax[1].set_xlim((0,is_match.shape[1]))
  794. ax[1].set_xlabel('Cell (#)')
  796. ax[2].scatter(consistency_iter*100,np.arange(num_iter),color='k',s=1,edgecolor='k',marker='.')
  797. ax[2].axvline(1/k*100,color='k',linestyle=':',linewidth=.25)
  798. ax[2].set_xlabel('Iteration\nConsistency (%)',va='bottom')
  799. ax[2].set_xlim((102,0))
  800. ax[2].set_ylim((0,num_iter))
  801. ax[2].yaxis.tick_right()
  802. ax[2].xaxis.set_label_position('top')
  803. ax[2].xaxis.tick_top()
  806. # fig.delaxes(ax[3])
  807. ax[3].errorbar(mean_consistency_iter*100, mean_consistency_cell*100, yerr=sd_consistency_cell*100, xerr=std_consistency_iter*100,color='k',capsize=.75,linewidth=0.5,capthick=0.5)
  808. ax[3].axvline(1/k*100,color='k',linestyle=':',linewidth=.25)
  809. ax[3].axhline(1/k*100,color='k',linestyle=':',linewidth=.25)
  810. ax[3].set_xlim((100,0))
  811. ax[3].set_ylim((0,100))
  812. ax[3].yaxis.tick_right()
  814. # fig_cv_prob.suptitle(f'k = {k} clusters', y=1.3,fontsize=plt.rcParams['font.size'])# y=1.075
  816. # fig_cv_prob.tight_layout(pad=.1)
  817. # fig_cv_prob.set_figwidth(1.8)
  818. # fig_cv_prob.set_figheight(1.6)
  819. ax[0].set_xticks([])
  821. fig_cv_prob.patch.set_facecolor('white')
  823. plt.show()
  825. plt.tight_layout()
  827. fig_cv_prob.savefig(f'{folder}/CrossVal_Data k={k} Labeling.'+fig_opt['format'],**fig_opt)
  831. return cval_res
  836. # %%
  837. ### Git clone data tables for sorting ###
  838. try: shutil.rmtree('./'+git_rep)
  839. except: None
  840. git_rep = 'Prelim_2023'
  841. git_link = 'https://github.com/dtabuena/'+git_rep+'/'
  842. !git clone $git_link
  843. file_list = [f for f in [f for r,d,f in os.walk('./'+git_rep+'/Cleaned_DataSets/' )][0] if '.xlsx' in f]
  844. print(file_list)
  845. mpl.font_manager.fontManager.addfont('./'+git_rep+'/Cleaned_DataSets/arial.ttf')
  847. # %%
  848. fig_cartoon_x(fig_opt={'dpi': 300, 'format':'svg','bbox_inches':None})
  850. # %%
  851. ### Cycle throught Data tables, cluster, and plot. ###
  855. multipack={}
  857. fig_opt = {'dpi': 300, 'format':'svg','bbox_inches':None}
  860. for f in file_list: # [f for f in file_list if 'CA3' in f]:#
  861. my_df_dict, xl_to_analyze = read_data_file(f,dir='./'+git_rep+'/Cleaned_DataSets/')
  862. folder = xl_to_analyze.split('.')[0].replace(' ','_')
  863. try: os.makedirs(folder)
  864. except: None
  865. full_df = dict_to_df(my_df_dict,AP_cut=60)
  866. labeling = ['cell','Cell','type','extra']
  867. numericals, data_for_fit, full_df_no_null = clean_data(labeling, full_df)
  868. type_num, num_t, uniq_types,type_list = get_types(full_df_no_null)
  869. packaged = cluster_full(data_for_fit,xl_to_analyze,numericals)
  870. packaged['type_num'] = type_num
  871. if 'CA1' in f: packaged['cell_type'] ='CA1 PC'
  872. if 'CA3' in f: packaged['cell_type'] ='CA3 PC'
  873. if 'II' in f: packaged['cell_type'] ='DGC II'
  874. if ' I ' in f: packaged['cell_type'] ='DGC I'
  875. df_for_34_2W_anova, tukey_data = do_stats(type_list,packaged['score'])
  876. plot_summary(num_t,numericals,type_list,uniq_types,packaged,df_for_34_2W_anova,tukey_data,folder,fig_opt =fig_opt)
  877. k2_cval = cross_val(data_for_fit,numericals,packaged['new_cent_order'],folder,k=2,num_iter = 1000,SubSampleFrac = .8,
  878. cell_type=packaged['cell_type'],fig_opt =fig_opt,ds_iter=1)
  879. k3_cval = cross_val(data_for_fit,numericals,packaged['new_cent_order'],folder,k=3,num_iter = 1000,SubSampleFrac = .8,
  880. cell_type=packaged['cell_type'],fig_opt =fig_opt,ds_iter=1)
  881. packaged['k2_cval'] = k2_cval
  882. packaged['k3_cval'] = k3_cval
  883. multipack[packaged['cell_type']] = packaged
  884. save_dl(folder)
  886. # %%
  889. # %%
  890. ########### Write K-Means Values ################
  891. all_dfs = dict()
  892. for celltype, data in multipack.items():
  893. # print(celltype)
  894. # print(data.keys())
  895. # print(data['features'])
  896. # print(data['k2_cval'])
  898. # Create dictionary with features column
  899. df_dict = {'features': data['features']}
  901. # Add a column for each cluster center
  902. for i in range(data['centers'].shape[1]):
  903. df_dict[f'full_centroid_K{i}'] = data['centers'][:, i]
  905. means = data['k2_cval']['mean_center']
  906. for i in range(means.shape[1]):
  907. df_dict[f'mean_centroid_K{i}'] = means[:, i]
  909. confidence = data['k2_cval']['confidence']
  910. for i in range(confidence.shape[2]): # Iterate over clusters (last dimension)
  911. df_dict[f'confidence_low_K{i}'] = confidence[0, :, i] # Low bound
  912. df_dict[f'confidence_high_K{i}'] = confidence[1, :, i] # High bound
  914. feat_df = pd.DataFrame(df_dict)
  915. feat_df = feat_df.sort_values('full_centroid_K0', ascending=False)
  916. feat_df = feat_df.rename(columns=lambda x: x.replace('K0', 'Hyper').replace('K1', 'Normal'))
  917. all_dfs[celltype] = feat_df
  920. # Write to multi-tabbed Excel file
  921. with pd.ExcelWriter('kmeans_results.xlsx', engine='openpyxl') as writer:
  922. for celltype, df in all_dfs.items():
  923. df.to_excel(writer, sheet_name=celltype, index=False)
  925. print("Excel file 'kmeans_results.xlsx' created successfully!")
  926. colab.files.download('kmeans_results.xlsx')
  928. # %%
  929. ############### Combined Centers Plot ####################
  930. FS=6
  931. # set_font_all(FS)
  932. dpi=300
  933. num_types = len(multipack)
  934. fig_combined_centers, ax = plt.subplots(1,1,figsize=np.array([7.25,1.6])*1,dpi=dpi,layout="constrained")
  935. cmap = plt.cm.cool(np.linspace(0,1,2))*np.array([0.6]*3+[1])
  936. dark_cool = mpl.colors.ListedColormap(cmap)
  937. # cent_ord = multipack['CA3 PC']['new_cent_order']
  938. features = np.array(multipack['CA3 PC']['features'])
  940. mean_center = multipack['CA3 PC']['k2_cval']['mean_center']
  941. plot_order = np.argsort(np.diff(mean_center,1).T)[0]
  943. keys = multipack.keys()
  944. keys = ['CA3 PC', 'DGC II', 'DGC I', 'CA1 PC' ]
  945. for ci in range(num_types):
  946. # ci_type = list(multipack.keys())[ci]
  947. ci_type = keys[ci]
  949. centers = multipack[ci_type]['centers']
  951. bar_wid = 1/(num_types*2+2)
  952. x = np.arange(len(features))
  954. m_color = cmap[0,:]*[1,1,1,(ci+1)/4]
  955. m_ys = centers[:,0]
  956. m_label = ci_type + ' - Hyper'
  957. m_x = x + (ci+1)*bar_wid - 0.5
  958. ax.bar(m_x,m_ys[plot_order],color=m_color,edgecolor='k',width=bar_wid,label=m_label,linewidth=.5)
  960. l_color = cmap[1,:]*[1,1,1,(ci+1)/4]
  961. l_ys = centers[:,1]
  962. l_label = ci_type + ' - Normal'
  963. l_x = m_x + (num_types)*bar_wid
  964. ax.bar(l_x,l_ys[plot_order],color=l_color,edgecolor='k',width=bar_wid,label=l_label,linewidth=.5)
  966. ax.set_xticks(x)
  967. ax.set_ylim(bottom=-1, top=1.3)
  968. ax.set_xlim(-0.5,x[-1]+.5)
  969. ax.set_xticklabels(features[plot_order],rotation=25,ha='right', rotation_mode='anchor')
  970. ax.axhline(0,color='k',linewidth=.5)
  971. # ax.axvline(0,color='k')
  972. ax.set_ylabel('z-scored'+' '+'values')
  973. # ax.title.set_text('Cluster Centers')
  974. ax.text(np.mean(ax.get_xlim()),ax.get_ylim()[-1]*.95,'Cluster Centers',ha='center',va='top')
  975. handles, bar_labels = plt.gca().get_legend_handles_labels()
  976. order = np.arange(8)
  977. order = np.append(np.arange(num_types)*2,np.arange(num_types)*2+1)
  978. print(order)
  979. ax.legend([handles[i] for i in order],[bar_labels[i] for i in order],ncol=1,loc='center left',frameon=False,bbox_to_anchor=(1, .5)) ## ,bbox_to_anchor=(1, 0), loc='upper center'
  980. # plt.()
  981. plt.show()
  982. fig_combined_centers.savefig('combined_centers.'+fig_opt['format'],**fig_opt)
  983. colab.files.download('combined_centers.'+fig_opt['format'])
  988. # %%
  990. from cmh import CMH
  992. print(type_list)
  996. ##### DG TypeR ######
  997. type_list =['DG Type II', 'DG Type I']
  998. cmh_df = pd.DataFrame(columns=['Geno','Age','Type'])
  999. counts = {}
  1001. e4_table =pd.DataFrame(columns=type_list,index=['7-9mo','17-19mo'])
  1002. e3_table = e4_table.copy()
  1003. for f in file_list:
  1004. if 'DG' in f:
  1005. c_type = [t for t in type_list if t+' ' in f][0]
  1006. my_df_dict,_ = read_data_file(f)
  1007. for k,v in my_df_dict.items():
  1008. if '19' in k: age = '17-19mo'
  1009. else: age = '7-9mo'
  1010. if 'E4' in k: e4_table.at[age,c_type]=v.shape[0]
  1011. else: e3_table.at[age,c_type]=v.shape[0]
  1013. a = np.array([[k.split(' ')[0], k.split(' ')[1], c_type]]*v.shape[0])
  1014. df2 = pd.DataFrame(data=a, columns=['Geno','Age','Type'])
  1015. # cmh_df = cmh_df.append(df2, ignore_index=True) # DEPRICATED
  1016. cmh_df = pd.concat([cmh_df,df2], ignore_index=True) # REPLACED 2024.06.13
  1018. print('E3')
  1019. display(e3_table)
  1021. print('E4')
  1022. display(e4_table)
  1024. result = CMH(cmh_df, 'Type', 'Age', stratifier='Geno')
  1025. display(result)
  1027. # %%
  1028. def plot_summary_DG_comb(num_t,numericals,type_list,uniq_types,packaged,df_for_34_2W_anova,tukey_data,folder, clust_names = ['Hyper', 'Normal'],fig_opt = {'dpi': 300, 'format':'png',
  1029. 'bbox_inches':None}):
  1030. # fig_opt = {'dpi': 300, 'format':'svg','bbox_inches':None}
  1032. cmap = plt.cm.cool(np.linspace(0,1,2))*np.array([0.6]*3+[1])
  1033. dark_cool = mpl.colors.ListedColormap(cmap)
  1036. kmeans = packaged['kmeans']
  1037. centers = packaged['centers']
  1038. distances = packaged['distances']
  1039. rel_centers = packaged['rel_centers']
  1040. rel_counts = packaged['rel_counts']
  1041. total_D = packaged['total_D']
  1042. closeness = packaged['closeness']
  1043. distances_norm = packaged['distances_norm']
  1044. score = packaged['score']
  1045. data_for_fit = packaged['data_for_fit']
  1046. cell_type = packaged['cell_type']
  1047. FS=rcParams['font.size']
  1050. colors = plt.cm.viridis(np.linspace(0,1,num_t))
  1051. alpha_1 = 1.0
  1052. alpha_2 = 0.6
  1053. colors = [[4/255, 51/255, 255/255, alpha_1], [4/255, 51/255, 255/255, alpha_2],
  1054. [255/255, 38/255, 0/255, alpha_1], [255/255, 38/255, 0/255, alpha_2],
  1055. [0, 0.75, 0.75, alpha_1], [0, 0.75, 0.75, alpha_2],
  1056. [1, 0, 0.65, alpha_1], [1, 0, 0.65, alpha_2],
  1057. ]
  1059. list_colors = [colors[i] for i in type_num ]
  1060. color_dc={'apoE3-KI Type_I 7-9mo': colors[0],#[:-1]+[alpha_1],
  1061. 'apoE3-KI Type_II 7-9mo': colors[4],#[:-1]+[alpha_1],
  1062. 'apoE4-KI Type_I 7-9mo': colors[1],#[:-1]+[alpha_2],
  1063. 'apoE4-KI Type_II 7-9mo': colors[5],#[:-1]+[alpha_2],
  1064. 'apoE3-KI Type_I 17-19mo': colors[2],#[:-1]+[alpha_1],
  1065. 'apoE3-KI Type_II 17-19mo': colors[6],#[:-1]+[alpha_1],
  1066. 'apoE4-KI Type_I 17-19mo': colors[3],#[:-1]+[alpha_2],
  1067. 'apoE4-KI Type_II 17-19mo': colors[7],#[:-1]+[alpha_2],
  1068. }
  1072. '''
  1073. apoE3-KI Type_I 7-9mo
  1074. apoE3-KI Type_II 7-9mo
  1075. apoE4-KI Type_I 7-9mo
  1076. apoE4-KI Type_II 7-9mo
  1077. apoE3-KI Type_I 17-19mo
  1078. apoE3-KI Type_II 17-19mo
  1079. apoE4-KI Type_I 17-19mo
  1080. apoE4-KI Type_II 17-19mo
  1081. '''
  1083. '''raw_data_map '''
  1084. counts = [np.sum( [tt==t for tt in type_list]) for t in uniq_types ]
  1085. fig_norm_vals, ax = plt.subplots(1,num_t,figsize = [4.3,1],gridspec_kw={'width_ratios': counts},layout='constrained')
  1086. data_X_min = np.min(data_for_fit)
  1087. data_X_max = np.max(data_for_fit)
  1088. for ti in range(num_t):
  1089. t = uniq_types[ti]
  1090. t_colr = colors[ti]
  1091. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  1092. c = ax[ti].pcolor(data_for_fit.T[:,is_type], vmin=data_X_min*1.2, vmax=data_X_max*1.2)
  1093. ax[ti].set_yticks(np.arange(len(numericals))+0.5)
  1094. if ti ==0 : ax[ti].set_yticklabels(numericals,rotation=0)
  1095. else: ax[ti].get_yaxis().set_visible(False)
  1096. ax[ti].set_xlabel('Cell (#)')
  1097. ax[ti].title.set_text(t.replace(' ','\n'))
  1098. start = min(is_type)+3 / len(is_type)
  1099. stop = max(is_type)-3 / len(is_type)
  1100. fig_norm_vals.suptitle(cell_type+ ' Normalized Parameters',y=1.05, fontsize=FS)
  1101. # plt.tight_layout()
  1102. plt.colorbar(c, ax=ax[ti],label='z-scored'+'\n'+'values')
  1103. plt.show()
  1106. ################# Centers
  1107. new_cent_order = np.argsort(centers[:,1]).tolist()[::-1]
  1108. # fig_clust_pie, ax = plt.subplots(1,2,figsize = [12,5], gridspec_kw={'width_ratios': [2, 1]} )
  1109. fig_centers, ax = plt.subplots(1,1,figsize = [3,2],layout='constrained')
  1110. ax.set_prop_cycle(color=[cmap[0], cmap[1]])
  1111. ax.plot(centers[new_cent_order,:],'-o',linewidth=.5)
  1112. ax.legend(clust_names,loc='lower center',frameon=False)
  1113. sorted_numericals = [numericals[i] for i in new_cent_order]
  1114. ax.set_xticks(np.arange(len(numericals)))
  1115. ax.set_xticklabels(sorted_numericals,rotation=45,ha='right', rotation_mode='anchor',fontsize=FS)
  1116. ax.set_ylabel('z-scored'+' '+'values')
  1117. ax.axhline(0,color='k')
  1118. ax.title.set_text(cell_type +' '+'Cluster Centers')
  1119. # plt.tight_layout()
  1121. ################# Split Distance heat map
  1122. counts = [np.sum( [tt==t for tt in type_list]) for t in uniq_types ]
  1123. print(counts)
  1124. b_dim = 2
  1125. # fig_distances_split, ax = plt.subplots(1,num_t,figsize = [8,3],gridspec_kw={'width_ratios': counts},constrained_layout=True) ## OG
  1126. fig_distances_split, ax = plt.subplots(1,num_t,figsize = [7.25,.5],gridspec_kw={'width_ratios': counts},layout='constrained') ## NEw
  1127. dist_min = np.min(distances_norm)
  1128. dist_max = np.max(distances_norm)
  1129. for ti in range(num_t):
  1130. t = uniq_types[ti]
  1131. t_colr = colors[ti]
  1132. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  1133. c = ax[ti].pcolor(distances_norm.T[:,is_type], vmin=dist_min*1.2, vmax=dist_max*1.2)
  1134. ax[ti].set_yticks(np.arange(len(clust_names))+0.5)
  1135. clust_names_nl = [c.replace(' ',' ') for c in clust_names]
  1136. if ti ==0 : ax[ti].set_yticklabels(clust_names_nl,rotation=0)
  1137. else: ax[ti].get_yaxis().set_visible(False)
  1138. # ax[ti].title.set_text(t.replace(' ','\n'))
  1139. # ax[ti].title.set_fontsize(16)
  1140. # ax[ti].title.set_rotation(60)
  1141. ax[ti].set_xlabel('Cell (#)')
  1142. # fig_distances_split.tight_layout()
  1143. fig_distances_split.colorbar(c, ax=ax[ti],label='z-scored'+'\n'+'distance')
  1144. for ti in range(num_t):
  1145. t = uniq_types[ti]
  1146. ax[ti].text(0.5, 1, t.replace(' ','\n'),transform=ax[ti].transAxes,ha='center',rotation=0,va='bottom',fontsize=FS)
  1147. plt.show()
  1148. # xxxxxxx
  1150. ############# Combined Scatter
  1151. fig_comb_scatter, ax = plt.subplots(1,figsize = [4,2],layout='constrained')
  1152. for ti in range(num_t):
  1153. t = uniq_types[ti]
  1154. t_colr = color_dc[t]
  1155. print(t,t_colr)
  1156. t_colr_edge = t_colr[:-1] + [alpha_1]
  1157. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  1158. ax.scatter(distances_norm[is_type,0],distances_norm[is_type,1],color=t_colr, label=t,linewidth=1,edgecolor=t_colr_edge,s=4)
  1159. ax.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
  1160. # ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
  1161. ax.set_xlabel('Dist. from '+clust_names[0]+' Center \n (z-scored)')
  1162. ax.set_ylabel('Dist. from '+clust_names[1]+' Center \n (z-scored)')
  1163. ax.plot([0, 1], [0, 1], transform=ax.transAxes,color='k')
  1164. ax.text(0.98, 0.02, ''+clust_names[1],transform=ax.transAxes,ha='right',va='bottom',fontsize=FS)
  1165. ax.text(0.02, 0.98, ''+clust_names[0],transform=ax.transAxes,ha='left',va='top',fontsize=FS)
  1166. comb_xlim = ax.get_xlim()
  1167. comb_ylim = ax.get_ylim()
  1168. # print(comb_xlim)
  1169. # print(comb_ylim)
  1170. common_lim = np.array([-1,1]) * np.max(np.abs(list(comb_xlim)+list(comb_ylim)))
  1171. common_lim= [np.min([comb_xlim[0],comb_ylim[0]]), np.max([comb_xlim[1],comb_ylim[1]]) ]
  1172. # print(common_lim)
  1173. ax.set_xlim(common_lim)
  1174. ax.set_ylim(common_lim)
  1176. #### separate Scatter #####
  1177. b_dim = 1.5
  1178. if 'CA3' in cell_type: fig_sep_scatter, axs = plt.subplots(1,num_t,figsize = [num_t*b_dim,b_dim*1.2],layout='constrained')
  1179. else: fig_sep_scatter, axs = plt.subplots(2,int(num_t/2),figsize = [int(num_t/2)*b_dim,b_dim*2*1.2],layout='constrained') # [int(num_t/2)*b_dim,b_dim*2]
  1180. axs = axs.T.flatten()
  1181. for ti in range(num_t):
  1182. t = uniq_types[ti]
  1183. t_colr = color_dc[t]
  1184. t_colr_edge = t_colr[:-1] + [alpha_1]
  1185. is_type = [ind for ind in range(len(type_num)) if type_num[ind]==ti]
  1186. axs[ti].scatter(distances_norm[is_type,0],distances_norm[is_type,1],color=t_colr, label=t,linewidth=1,edgecolor=t_colr_edge,s=36)
  1187. axs[ti].plot([0, 1], [0, 1], transform=axs[ti].transAxes,color=[0,0,0,1])
  1188. axs[ti].title.set_text(t.replace(' ','\n'))
  1189. axs[ti].set_xlim(common_lim)
  1190. axs[ti].set_ylim(common_lim)
  1193. axs[ti].set_xlabel(clust_names[0]+' Dist.',fontsize=FS)
  1194. axs[ti].set_ylabel(clust_names[1]+' Dist.',fontsize=FS)
  1195. axs[ti].text(0.98, 0.02, ''+clust_names[1],transform=axs[ti].transAxes,ha='right',va='bottom',fontsize=FS)
  1196. axs[ti].text(0.02, 0.98, ''+clust_names[0],transform=axs[ti].transAxes,ha='left',va='top',fontsize=FS)
  1198. x=distances_norm[is_type,0]
  1199. y=distances_norm[is_type,1]
  1200. for i in range(len(x)):
  1201. plt.plot([x[i], x[i]], [y[i], x[i]], 'k:')
  1203. # plt.tight_layout()
  1206. ######## ''' Violins of Residual Plots '''
  1207. fig_violin_residuals, ax = plt.subplots(1,figsize=[7.25*.5,1.4],layout='constrained')
  1208. comb_list = list(set(df_for_34_2W_anova['combination']))
  1210. young_types = [t for t in comb_list if '-9' in t and 'Cre' not in t]
  1211. young_types.sort()
  1212. cre_types = [t for t in comb_list if 'Cre' in t]
  1213. cre_types.sort()
  1214. old_types = [t for t in comb_list if '-19' in t and 'Cre' not in t]
  1215. old_types.sort()
  1216. comb_list = young_types + cre_types + old_types
  1218. sns_dict={}
  1219. print('comb_list',comb_list)
  1220. for ti in range(len(comb_list)):
  1221. type_i = comb_list[ti]
  1222. match_type = df_for_34_2W_anova['combination']==type_i
  1223. y_vals = list(np.array(df_for_34_2W_anova[match_type]['T1_Score']))
  1224. x_vals = np.ones_like(y_vals)*ti + int(np.floor(ti/2))
  1225. pos = list(np.ones_like(y_vals)*ti)
  1226. t_colr = colors[ti]
  1227. t_colr = color_dc[type_i]
  1228. sns_dict[type_i] =y_vals
  1230. print(sns_dict.keys())
  1231. sns_df = pd.DataFrame.from_dict(sns_dict, orient='index').T
  1232. display(sns_df.head())
  1234. colors_sns = [color_dc[t] for t in sns_df.columns]
  1235. colors_sns = [ list(np.array(c[:-1]) + (1 - np.array(c[:-1]))* (1-c[-1])) + list(c[-1:]) for c in colors_sns]
  1237. print("colors_sns",colors_sns,"colors_sns")
  1239. # Plot the swarm plot
  1240. sns.swarmplot(ax=ax, data=sns_df, palette=colors_sns, label=t, linewidth=.05, edgecolor='k', size=2)
  1242. # Plot the violin plot
  1243. sns.violinplot(ax=ax, data=sns_df, palette=colors_sns, scale='area', label=t, inner=None, linewidth=0)
  1245. # Adjust the violins to set the face and edge colors
  1246. violins = [v for v in ax.collections if isinstance(v, matplotlib.collections.PolyCollection)] # Get all violins
  1248. # Iterate over violins and apply fill and edge color
  1249. for i, violin in enumerate(violins):
  1250. violin.set_facecolor([0.9, 0.9, 0.9]) # Set fill color to light gray
  1251. violin.set_edgecolor(colors_sns[i % len(colors_sns)]) # Set edge color from colors_sns
  1252. violin.set_linewidth(1) # Set the linewidth to ensure the edges are visible
  1254. ax.text(0, 1, clust_names[0],transform=ax.transAxes,ha='right',va='center',fontsize = FS,color=[.6, 0, .6])
  1255. ax.text(0, 0, clust_names[1],transform=ax.transAxes,ha='right',va='center',fontsize = FS,color=[0, .6, .6]) # =
  1257. ax.set_ylabel('Residual Value')
  1258. ax.set_ylim(-7.5,9.5)
  1259. ax.axhline(0,color='k')
  1260. ax_labels = sns_df.columns
  1261. ax.set_xticks(ticks=range(len(ax_labels)))
  1262. ax_labels_r = [c.replace(' ','\n') for c in ax_labels]
  1263. ax_labels_r = [c.replace('\nsyn',' syn') for c in ax_labels_r]
  1264. ax_labels_r = [c.replace('\n7-9mo','') for c in ax_labels_r]
  1265. ax_labels_r = [c.replace('\n17-19mo','') for c in ax_labels_r]
  1266. ax_labels_r = [c.replace('-Cre', '\nCre') for c in ax_labels_r]
  1267. ax_labels_r = [c.replace('_', ' ') for c in ax_labels_r]
  1268. ax_labels_r = [c.replace('apo', '') for c in ax_labels_r]
  1269. print('ax_labels',ax_labels)
  1270. print('ax_labels_r',ax_labels_r)
  1271. ax.set_xticklabels(ax_labels_r,rotation=0,ha='center',va='top', rotation_mode='anchor',fontsize=FS)
  1272. ax.title.set_text(cell_type+' Cluster Membership')
  1274. div_x = np.sum( ['7-9' in l for l in ax_labels])-0.5
  1275. segment_x = [div_x, div_x]
  1276. # segment_y = ax.get_ylim()*np.array([1.8,1])
  1277. segment_y = ax.get_ylim()*np.array([1,1])
  1278. print('segments',segment_x,segment_y)
  1279. ll = plt.plot(segment_x,segment_y,':k') # returns a sequence of line objects
  1280. ll[0].set_clip_on(False)
  1282. month_offset = 0
  1283. plt.text(div_x/2-0.25, segment_y[0]-month_offset, '7-9mo', fontsize=FS,ha='center',va='bottom')
  1284. plt.text(div_x+div_x/2+0.25, segment_y[0]-month_offset, '17-19mo', fontsize=FS,ha='center',va='bottom')
  1287. labels = comb_list
  1288. true_tuckey = tukey_data[tukey_data['reject']==True]
  1290. astk_range = [0.05,0.01,0.001,0.0001]
  1291. for r in true_tuckey.index:
  1292. g1 = true_tuckey.loc[r,'group1']
  1293. g2 = true_tuckey.loc[r,'group2']
  1294. g1 = g1.replace('_II','X')
  1295. g1 = g1.replace('_I','X')
  1296. g2 = g2.replace('_II','X')
  1297. g2 = g2.replace('_I','X')
  1298. is_match = g1 in g2
  1301. if is_match:
  1302. start_l = true_tuckey.loc[r,'group1']
  1303. stop_l = true_tuckey.loc[r,'group2']
  1304. start_x = labels.index(start_l)
  1305. stop_x = labels.index(stop_l)
  1306. y0 = np.max(df_for_34_2W_anova['T1_Score'])*1.25
  1307. y_offset = abs((stop_x-start_x))
  1308. y = y0+ y_offset * .75 + np.min([stop_x,start_x])*.25
  1310. half_dx = abs(stop_x-start_x)/2
  1311. ax.errorbar( np.mean([stop_x,start_x]), y,yerr=0, xerr=half_dx,color='k',capsize=.75,linewidth=0.75,capthick=0.75)
  1313. cent =np.mean([start_x,stop_x])
  1314. num_ask = np.sum([true_tuckey.loc[r,'p-adj'] < t for t in astk_range] )
  1316. plt.text(cent,y, '*' * num_ask,fontsize=FS+10,ha='center',va='center')
  1319. fig_size = fig_violin_residuals.get_size_inches() # Get the figure size
  1321. fig_norm_vals.savefig(f'{folder}/Normalized_Vals_split.'+fig_opt['format'],**fig_opt)
  1322. fig_violin_residuals.savefig(f'{folder}/fig_violin_residuals.'+fig_opt['format'],**fig_opt)
  1323. fig_distances_split.savefig(f'{folder}/fig_distances_split.'+fig_opt['format'],**fig_opt)
  1324. fig_sep_scatter.savefig(f'{folder}/GenoType_dist_Separate.'+fig_opt['format'],**fig_opt)
  1325. fig_comb_scatter.savefig(f'{folder}/GenoType_dist.'+fig_opt['format'],**fig_opt)
  1326. return None
  1329. ### testing
  1330. # plot_summary_DG_comb(num_t,numericals,type_list,uniq_types,packaged,df_for_34_2W_anova,tukey_data,folder,fig_opt =fig_opt)
  1332. # %%
  1333. ##### DG Combined ####
  1334. file_list_DG = [f for f in file_list if 'DG' in f]
  1335. print(file_list_DG)
  1336. dg_dict={}
  1337. for f in file_list_DG:
  1338. f_dict,_=read_data_file(f)
  1339. for k,v in f_dict.items():
  1340. if 'II' in f: k=k.replace('KI','KI Type_II')
  1341. if ' I ' in f: k=k.replace('KI','KI Type_I')
  1342. dg_dict[k]=v
  1344. folder ='DG_Comb'
  1345. try: os.makedirs(folder)
  1346. except: None
  1347. full_df = dict_to_df(dg_dict,AP_cut=60)
  1348. labeling = ['cell','Cell','type','extra']
  1349. numericals, data_for_fit, full_df_no_null = clean_data(labeling, full_df)
  1350. type_num, num_t, uniq_types,type_list = get_types(full_df_no_null)
  1351. packaged = cluster_full(data_for_fit,folder,numericals)
  1352. packaged['cell_type']='Combined DGC'
  1353. df_for_34_2W_anova, tukey_data = do_stats(type_list,packaged['score'])
  1356. _ = [print(t)for t in uniq_types]
  1358. plot_summary_DG_comb(num_t,numericals,type_list,uniq_types,packaged,df_for_34_2W_anova,tukey_data,folder,fig_opt =fig_opt)
  1359. packaged['k2_cval'] = cross_val(data_for_fit,numericals,packaged['new_cent_order'],folder,k=2,num_iter = 1000,SubSampleFrac = .8,
  1360. cell_type=packaged['cell_type'],fig_opt =fig_opt,ds_iter=5)
  1361. # _ = cross_val(data_for_fit,numericals,packaged['new_cent_order'],folder,k=3,num_iter = 1000,SubSampleFrac = .8,
  1362. # cell_type=packaged['cell_type'],fig_opt =fig_opt,ds_iter=5)
  1363. # _ = cross_val(data_for_fit,numericals,packaged['new_cent_order'],folder,k=4,num_iter = 1000,SubSampleFrac = .8,
  1364. # cell_type=packaged['cell_type'],fig_opt =fig_opt,ds_iter=5)
  1365. # _ = cross_val(data_for_fit,numericals,packaged['new_cent_order'],folder,k=5,num_iter = 1000,SubSampleFrac = .8,
  1366. # cell_type=packaged['cell_type'],fig_opt =fig_opt,ds_iter=5)
  1367. save_dl(folder)
  1369. # %%
  1370. data=packaged
  1372. # Create dictionary with features column
  1373. df_dict = {'features': data['features']}
  1375. # Add a column for each cluster center
  1376. for i in range(data['centers'].shape[1]):
  1377. df_dict[f'full_centroid_K{i}'] = data['centers'][:, i]
  1379. means = data['k2_cval']['mean_center']
  1380. for i in range(means.shape[1]):
  1381. df_dict[f'mean_centroid_K{i}'] = means[:, i]
  1383. confidence = data['k2_cval']['confidence']
  1384. for i in range(confidence.shape[2]): # Iterate over clusters (last dimension)
  1385. df_dict[f'confidence_low_K{i}'] = confidence[0, :, i] # Low bound
  1386. df_dict[f'confidence_high_K{i}'] = confidence[1, :, i] # High bound
  1388. feat_df = pd.DataFrame(df_dict)
  1389. feat_df = feat_df.sort_values('full_centroid_K0', ascending=False)
  1390. feat_df = feat_df.rename(columns=lambda x: x.replace('K0', 'Hyper').replace('K1', 'Normal'))
  1391. feat_df = feat_df.set_index('features')
  1392. display(feat_df)
  1393. feat_df.to_csv('combined_dgcs.csv')
  1394. colab.files.download('combined_dgcs.csv')
  1396. # %%
  1397. ## NatAging Resubmission
  1398. """
  1399. New Supplemental Figure
  1400. """
  1402. row_0 =0
  1403. row_1 =row_0+160
  1404. row_2 =row_1+170
  1405. row_3 = row_2+190
  1406. row_4 = row_3+215
  1407. dpi = 96
  1408. width = 7.25*dpi
  1409. length = 9.4*dpi
  1410. scale_factor= 1.33
  1411. weight='bold'
  1412. label_size = 12
  1415. sup_layout_1 = Figure(str(width), str(length),
  1416. Panel(SVG("/content/CA3_parameters_for_PCA_Clean/CrossVal_Data k=2 Centers.svg"),
  1417. Text("a", 10, 20, size=label_size,weight=weight,font='arial')
  1418. ).scale(scale_factor).move(width*0, row_0),
  1419. Panel(SVG("/content/DG_Type_II_GC_main_parameters_Clean/CrossVal_Data k=2 Centers.svg"),
  1420. Text("b", 10, 20, size=label_size,weight=weight)
  1421. ).scale(scale_factor).move(width*.25, row_0),
  1422. Panel(SVG("/content/DG_Type_I_GC_main_parameters_Clean/CrossVal_Data k=2 Centers.svg"),
  1423. Text("c", 10, 20, size=label_size,weight=weight)
  1424. ).scale(scale_factor).move(width*.50, row_0),
  1425. Panel(SVG("/content/CA1_main_parameters_Clean/CrossVal_Data k=2 Centers.svg"),
  1426. Text("d", 10, 20, size=label_size,weight=weight)
  1427. ).scale(scale_factor).move(width*.75, row_0),
  1430. Panel(SVG("/content/CA3_parameters_for_PCA_Clean/CrossVal_Data k=2 Labeling.svg"),
  1431. Text("e", 10, 20, size=label_size, weight=weight)
  1432. ).scale(scale_factor).move(width*0, row_1),
  1433. Panel(SVG("/content/DG_Type_II_GC_main_parameters_Clean/CrossVal_Data k=2 Labeling.svg"),
  1434. Text("f", 10, 20, size=label_size, weight=weight)
  1435. ).scale(scale_factor).move(width*.25, row_1),
  1436. Panel(SVG("/content/DG_Type_I_GC_main_parameters_Clean/CrossVal_Data k=2 Labeling.svg"),
  1437. Text("g", 10, 20, size=label_size, weight=weight)
  1438. ).scale(scale_factor).move(width*.50, row_1),
  1439. Panel(SVG("/content/CA1_main_parameters_Clean/CrossVal_Data k=2 Labeling.svg"),
  1440. Text("h", 10, 20, size=label_size, weight=weight)
  1441. ).scale(scale_factor).move(width*.75, row_1),
  1443. Panel(SVG("/content/CA3_parameters_for_PCA_Clean/CrossVal_Data k=3 Labeling.svg"),
  1444. Text("i", 10, 20, size=label_size, weight=weight)
  1445. ).scale(scale_factor).move(width*.0, row_2),
  1446. Panel(SVG("/content/DG_Type_II_GC_main_parameters_Clean/CrossVal_Data k=3 Labeling.svg"),
  1447. Text("j", 10, 20, size=label_size, weight=weight)
  1448. ).scale(scale_factor).move(width*.25, row_2),
  1449. Panel(SVG("/content/DG_Type_I_GC_main_parameters_Clean/CrossVal_Data k=3 Labeling.svg"),
  1450. Text("k", 10, 20, size=label_size, weight=weight)
  1451. ).scale(scale_factor).move(width*.50, row_2),
  1452. Panel(SVG("/content/CA1_main_parameters_Clean/CrossVal_Data k=3 Labeling.svg"),
  1453. Text("l", 10, 20, size=label_size, weight=weight)
  1454. ).scale(scale_factor).move(width*.75, row_2),
  1456. Panel(SVG("/content/CA3_parameters_for_PCA_Clean/GenoType_dist_Separate.svg"),
  1457. Text("m", 10, 20, size=label_size, weight=weight)
  1458. ).scale(scale_factor).move(width*0-0, row_3),
  1459. Panel(SVG("/content/DG_Type_II_GC_main_parameters_Clean/GenoType_dist_Separate.svg"),
  1460. Text("n", 10, 20, size=label_size, weight=weight)
  1461. ).scale(scale_factor).move(width*.25, row_3),
  1462. Panel(SVG("/content/DG_Type_I_GC_main_parameters_Clean/GenoType_dist_Separate.svg"),
  1463. Text("o", 10, 20, size=label_size, weight=weight)
  1464. ).scale(scale_factor).move(width*.5, row_3),
  1465. Panel(SVG("/content/CA1_main_parameters_Clean/GenoType_dist_Separate.svg"),
  1466. Text("p", 10, 20, size=label_size, weight=weight)
  1467. ).scale(scale_factor).move(width*.75, row_3),
  1469. Panel(SVG("//content/DG_Comb/CrossVal_Data k=2 Centers.svg"),
  1470. Text("q", 10, 0, size=label_size,weight=weight)
  1471. ).scale(scale_factor).move(width*.25, row_4),
  1473. Panel(SVG("/content/DG_Comb/fig_violin_residuals.svg"),
  1474. Text("r", 10, 0, size=label_size,weight=weight)
  1475. ).scale(scale_factor).move(width*.495, row_4),
  1478. # Add vertical dividing lines
  1479. Line(points=[(width * 0.25+4.5, 15), (width * 0.25+4.5, length*.9)], width=1, color='grey'),
  1480. Line(points=[(width * 0.50+4.5, 15), (width * 0.5+4.5, length*.8)], width=1, color='grey'),
  1481. Line(points=[(width * 0.75+4.5, 15), (width * 0.75+4.5, length*.8)], width=1, color='grey'),
  1483. )
  1486. sup_layout_1.save(f"SupFig4_v2.svg")
  1487. display(sup_layout_1)
  1488. colab.files.download(f"SupFig4_v2.svg")
  1490. # %%
  1493. ######### Alt fig 5
  1495. left_margin = 0
  1496. row_0 = 0
  1497. row_0_5 =row_0+20
  1498. row_1 =row_0+140
  1499. row_2 =row_1+180
  1500. dpi = 96
  1501. width = 7.25*dpi
  1502. length = 9.25*dpi
  1503. scale_factor= 1.333
  1504. weight='bold'
  1505. label_size = 12
  1507. weight = 'bold'
  1508. # label_size = label_size*scale
  1509. #### from svgutils.compose import *
  1511. CONFIG = {
  1512. "svg.file_path": ".",
  1513. "figure.save_path": ".",
  1514. "image.file_path": ".",
  1515. "text.position": (0, 0),
  1516. "text.weight": "bold",
  1517. "text.font": "Arial",
  1518. }
  1521. layout_1 = Figure(width,length,
  1522. Panel(SVG("/content/fig_cartoon_x.svg"),
  1523. Text("a", 0, 12-20, size=label_size, weight=weight)
  1524. ).move(width*.0, row_0_5).scale(scale_factor),
  1526. Panel(SVG("/content/CA1_main_parameters_Clean/GenoType_dist.svg"),
  1527. Text("b", -10, 12, size=label_size, weight=weight)
  1528. ).scale(scale_factor).move(width*.78, row_0),
  1530. Panel(SVG("/content/CA3_parameters_for_PCA_Clean/fig_violin_residuals.svg"),
  1531. Text("c", 0, 10, size=label_size, weight=weight)
  1532. ).scale(scale_factor).move(width*0, row_1),
  1533. Panel(SVG("/content/DG_Type_II_GC_main_parameters_Clean/fig_violin_residuals.svg"),
  1534. Text("d", 0, 10, size=label_size, weight=weight)
  1535. ).scale(scale_factor).move(width*.3+width*(.7/3*0), row_1),
  1536. Panel(SVG("/content/DG_Type_I_GC_main_parameters_Clean/fig_violin_residuals.svg"),
  1537. Text("e", 0, 10, size=label_size, weight=weight)
  1538. ).scale(scale_factor).move(width*.3+width*(.7/3*1), row_1),
  1539. Panel(SVG("/content/CA1_main_parameters_Clean/fig_violin_residuals.svg"),
  1540. Text("f", 0, 10, size=label_size, weight=weight)
  1541. ).scale(scale_factor).move(width*.3+width*(.7/3*2), row_1),
  1543. Panel(SVG("/content/combined_centers.svg"),
  1544. Text("g", 0, 0, size=label_size, weight=weight)
  1545. ).scale(scale_factor).move(left_margin+0, row_2),
  1547. )
  1551. display(layout_1)
  1552. layout_1.save("Figure_5_v2.svg")
  1553. colab.files.download('Figure_5_v2.svg')

_Clustering.ipynb at commit 3fa24de, no license · at the source

Overview

Authors: Dennis R. Tabuena1, Sung-Soo Jang1, Brian Grone1,2, Oscar Yip1,3, Emily A. Aery Jones1, Jessica Blumenfeld1,4, Zherui Liang1,4, Rajkamalpreet S. Mann1, Yaqiao Li1,2, Deanna Necula1,4, Nicole Koutsodendris1,5, Antara Rao1,5, Leonardo Ding1, Alex R. Zhang1, Yanxia Hao1,2, Qin Xu1,2, Seo Yeon Yoon1, Samuel De Leon1, Yadong Huang1,2,3,4,5,6,7, Misha Zilberter1
ORCID iDs: Dennis R. Tabuena, Oscar Yip, Emily A. Aery Jones, Jessica Blumenfeld, Yaqiao Li, Nicole Koutsodendris, Antara Rao, Leonardo Ding, Yadong Huang, Misha Zilberter
  1. Gladstone Institute of Neurological Disease, Gladstone Institutes,San Francisco, CA USA
  2. Gladstone Center for Translational Advancement, Gladstone Institutes,San Francisco, CA USA
  3. Biomedical Sciences Program, University of California, San Francisco,San Francisco, CA USA
  4. Neuroscience Program, University of California, San Francisco,San Francisco, CA USA
  5. Developmental and Stem Cell Biology Program, University of California, San Francisco,San Francisco, CA USA
  6. Department of Neurology, University of California, San Francisco,San Francisco, CA USA
  7. Department of Pathology, University of California, San Francisco,San Francisco, CA USA
Institutions: Gladstone Institutes (United States); Biomedical Sciences Program, University of California, San Francisco,San Francisco, CA USA; University of California, San Francisco (United States)
Journal: Nature aging, volume 6, issue 4, pages 886-904
Dates: received 10 July 2024; accepted 18 February 2026; published online 3 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s43587-026-01096-0 · PMID 41933197 · PMCID PMC13099648 · OpenAlex W7148741231
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), Alzheimer's / dementia (population)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging
Keywords: Alzheimer's disease
MeSH: Alzheimer Disease*, Apolipoprotein E4*, Hippocampus*, Nerve Net*, Neurons*, Aging, Animals, Dentate Gyrus, Disease Models, Animal, Male, Mice, Mice, Transgenic (* major topic)
Topic: Alzheimer's disease research and treatments (Physiology, Medicine), according to OpenAlex
Funding: National Institute on Aging (R01AG061150, R01AG087323, R01AG092390, F32AG085961, R01AG055682, R01AG071697, P01AG073082)
Citations: cited by 7 papers (Europe PMC); 110 references in the paper
Research resources: Primary antibodies included anti-NeuN RRID:AB_11205592, RRID:AB_2340470, 1:100) and anti-mCherry RRID:AB_2571870, RRID:AB_2782993

Abstract

The full impact of APOE4 (apolipoprotein E4), the strongest genetic risk factor for Alzheimer’s disease (AD), on neuronal and network function remains unclear, particularly during early preclinical stages of disease. Here we show that young APOE4 knockin (E4-KI) mice exhibit hippocampal region-specific network hyperexcitability that predicts later cognitive deficits. This early phenotype arises from cell-type-specific subpopulations of smaller, hyperexcitable neurons and is eliminated by selective removal of neuronal APOE4. With aging, E4-KI mice develop granule cell hyperexcitability, progressive inhibitory dysfunction and excitation–inhibition imbalance in the dentate gyrus. Single-nucleus RNA sequencing with multilevel gene filtering reveals age-dependent and cell-type-specific transcriptional changes and identifies candidate mediators of early neuronal hyperexcitability, including Nell2. Targeted CRISPR interference knockdown of Nell2 rescues abnormal excitability, implicating Nell2 as a contributor to APOE4-driven dysfunction. Together, these findings define molecular and circuit mechanisms linking neuronal APOE4-induced early network impairment to AD pathogenesis with aging.

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

Repository

Its files are read in the Code ↔ Paper reader above.

dtabuena/TabuenaJangGrone

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 3fa24de246b4d0adbc772091812095e48b2c3929, 17 December 2025
Languages: Jupyter (1)
Size: 5 files, 1 script
Software Heritage: not archived
Found in: “Code availability”
Holds: 1 notebook
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), Pingouin (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
1 file
At the source: github.com/dtabuena/TabuenaJangGrone

Code availability

All code generated during this study is accessible upon reasonable request to the corresponding authors. k-means clustering code is available as a GitHub repository (https://github.com/dtabuena/TabuenaJangGrone).

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

Tracing map

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

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • neither the text of the paper nor the code itself.

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

Data

Datasets cited

Data availability

All data associated with this study and the information of used materials are available in the main text, Methods or Supplementary Information. The snRNA-seq datasets of E3-KI and E4-KI mice at different ages are used from our previous publication61 (GEO accession ID: GSE167497 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE167497)). The snRNA-seq datasets of fE-KI/Syn1-Cre mice generated during the study are available at the GEO (accession ID: GSE279550 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE279550)). Source data are provided with this paper.

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

Versions

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

Version 2, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 20 authors, 1 keyword, 12 MeSH terms, 1 funder, 108 references, 4 RRIDs.

Cite

This paper

Tabuena, D. R., Jang, S.-S., Grone, B., Yip, O., Aery Jones, E. A., Blumenfeld, J., Liang, Z., Mann, R. S., Li, Y., Necula, D., Koutsodendris, N., Rao, A., Ding, L., Zhang, A. R., Hao, Y., Xu, Q., Yoon, S. Y., De Leon, S., Huang, Y., & Zilberter, M. (2026). Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis. Nature aging, 6(4), 886-904. https://doi.org/10.1038/s43587-026-01096-0

BibTeX

@article{tabuena2026neuronal,
author = {Tabuena, Dennis R. and Jang, Sung-Soo and Grone, Brian and Yip, Oscar and Aery Jones, Emily A. and Blumenfeld, Jessica and Liang, Zherui and Mann, Rajkamalpreet S. and Li, Yaqiao and Necula, Deanna and Koutsodendris, Nicole and Rao, Antara and Ding, Leonardo and Zhang, Alex R. and Hao, Yanxia and Xu, Qin and Yoon, Seo Yeon and De Leon, Samuel and Huang, Yadong and Zilberter, Misha},
title = {{Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis}},
journal = {Nature aging},
year = {2026},
month = apr,
volume = {6},
number = {4},
pages = {886--904},
publisher = {Nature Portfolio},
issn = {2662-8465},
doi = {10.1038/s43587-026-01096-0},
url = {https://doi.org/10.1038/s43587-026-01096-0},
pmid = {41933197},
pmcid = {PMC13099648}
}

RIS

TY - JOUR
AU - Tabuena, Dennis R.
AU - Jang, Sung-Soo
AU - Grone, Brian
AU - Yip, Oscar
AU - Aery Jones, Emily A.
AU - Blumenfeld, Jessica
AU - Liang, Zherui
AU - Mann, Rajkamalpreet S.
AU - Li, Yaqiao
AU - Necula, Deanna
AU - Koutsodendris, Nicole
AU - Rao, Antara
AU - Ding, Leonardo
AU - Zhang, Alex R.
AU - Hao, Yanxia
AU - Xu, Qin
AU - Yoon, Seo Yeon
AU - De Leon, Samuel
AU - Huang, Yadong
AU - Zilberter, Misha
TI - Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis
T2 - Nature aging
J2 - Nat Aging
PY - 2026
DA - 2026/04/03
VL - 6
IS - 4
SP - 886
EP - 904
SN - 2662-8465
PB - Nature Portfolio
DO - 10.1038/s43587-026-01096-0
UR - https://doi.org/10.1038/s43587-026-01096-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s43587-026-01096-0",
"type": "article-journal",
"title": "Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis",
"container-title": "Nature aging",
"author": [
{
"family": "Tabuena",
"given": "Dennis R."
},
{
"family": "Jang",
"given": "Sung-Soo"
},
{
"family": "Grone",
"given": "Brian"
},
{
"family": "Yip",
"given": "Oscar"
},
{
"family": "Aery Jones",
"given": "Emily A."
},
{
"family": "Blumenfeld",
"given": "Jessica"
},
{
"family": "Liang",
"given": "Zherui"
},
{
"family": "Mann",
"given": "Rajkamalpreet S."
},
{
"family": "Li",
"given": "Yaqiao"
},
{
"family": "Necula",
"given": "Deanna"
},
{
"family": "Koutsodendris",
"given": "Nicole"
},
{
"family": "Rao",
"given": "Antara"
},
{
"family": "Ding",
"given": "Leonardo"
},
{
"family": "Zhang",
"given": "Alex R."
},
{
"family": "Hao",
"given": "Yanxia"
},
{
"family": "Xu",
"given": "Qin"
},
{
"family": "Yoon",
"given": "Seo Yeon"
},
{
"family": "De Leon",
"given": "Samuel"
},
{
"family": "Huang",
"given": "Yadong"
},
{
"family": "Zilberter",
"given": "Misha"
}
],
"container-title-short": "Nat Aging",
"volume": "6",
"issue": "4",
"page": "886-904",
"DOI": "10.1038/s43587-026-01096-0",
"PMID": "41933197",
"PMCID": "PMC13099648",
"ISSN": "2662-8465",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s43587-026-01096-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
3
]
]
}
}

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.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: pandas, Matplotlib, NumPy, Alzheimer's / dementia, mouse, 10 references
[2] doi:10.1038/s41467-026-69866-3 [code]
Selective weakening of population-coupled synaptic activity in vivo in a mouse model of amyloid-beta pathology.
Journal: Nature communications
In common: Alzheimer's / dementia, mouse, 10 references
[3] doi:10.1016/j.celrep.2026.117646 [code]
Medial entorhinal-hippocampal desynchronization parallels the emergence of memory impairment in a mouse model of Alzheimer's disease pathology.
Journal: Cell reports
In common: pandas, SciPy, Matplotlib, 1 other tool, Alzheimer's / dementia, mouse, 7 references
[4] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: Pingouin, statsmodels, seaborn, 5 other tools, mouse, author Emily A. Aery Jones
[5] doi:10.1111/ejn.70480 [code]
Astrocyte Proximity Protects Synapses From Human Amyloid-Beta Induced Degeneration in a Mouse Ex Vivo Model of Early Alzheimer's Disease.
Journal: The European journal of neuroscience
In common: statsmodels, seaborn, scikit-learn, 4 other tools, Alzheimer's / dementia, mouse, 2 references
[6] doi:10.1038/s41467-026-74227-1 [code]
Age-related changes in behavioural and neural variability in a decision-making task.
Journal: Nature communications
In common: Pingouin, statsmodels, seaborn, 5 other tools, mouse, 1 reference
[7] doi:10.1111/acel.70494
Exceptional Longevity Modifying Allele APOE2 Promotes DNA Signaling Pathways Resisting Cellular Senescence in Human Neurons.
Journal: Aging cell
In common: Alzheimer's / dementia, 5 references
[8] doi:10.1371/journal.pcbi.1014571 [code]
SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: Pingouin, statsmodels, seaborn, 5 other tools, 1 reference
[9] doi:10.1126/sciadv.adz6517 [code]
Corticosterone-linked microglial activity underpins sexually dimorphic neuroplasticity after ketamine anesthesia.
Journal: Science advances
In common: statsmodels, seaborn, scikit-learn, 4 other tools, mouse, 2 references
[10] doi:10.3389/fnagi.2026.1847611 [code]
APOE ε4-associated hippocampal atrophy trajectories across the Alzheimer's disease continuum: a systematic review, meta-analysis, and longitudinal validation.
Journal: Frontiers in aging neuroscience
In common: statsmodels, seaborn, pandas, 3 other tools, Alzheimer's / dementia, 2 references

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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