Validation of human kallikrein 6 in the cerebrospinal fluid of patients with progressive and non-progressive alzheimer's disease: correlation with other biomarkers.

Code ↔ Paper

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

[1] § Results › Correlation of CSF hK6 concentration and APOE status ↔ KLK6_ELISA_data_analysis.ipynb, lines 189–199 · score 0.72 · e2e3, e2e4, e3e3, e3e4, status, APOE
[2] § Methods › Statistical analysis ↔ KLK6_ELISA_data_analysis.ipynb, lines 2178–2218 · score 0.64 · post hoc, Mann Whitney, Kruskal Wallis, Holm

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

Jupyter notebook · 2,452 lines · 98 KB · MIT · 2 matches

# %% [markdown]
# # Analysis of Kallikrein-6 protein (hK6) concentration in patients with Alzheimer's disease and controls
# %% [markdown]
# Important notes:
# Please read the README.md file prior to running the code to ensure reproducibility.
#
# The analysis herein is published in the following paper:
# Chatanaka MK, Soosaipillai A, Morato X, Diamandis E. Human Kallikrein 6 in the cerebrospinal fluid of patients with progressive and non-progressive Alzheimer’s disease. 2025. In preparation.
#
# %%
#important libraries to load
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.stats as stats # for stats.probplot, check QQ plot for normally distributed data
import seaborn as sns
from scipy.stats import f_oneway # one-way ANOVA
import statsmodels.api as sm
import statsmodels.formula.api as smf
import pingouin as pg #for for prtial correlation to adjust for sex
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multitest import multipletests #for post-hoc Holm correction
from dateutil.relativedelta import relativedelta #use this to accurately account for month and day differences when calculating age
# %%
#Read files based on location (SHOULD BE CHANGED BY USER)
file_path = "C:\\Users\\miyoh\\Git_Hub_Work\\Python_KLK6_data_analysis\\KLK6_ELISA_data_correct_for_dilution.xlsx" #Edit accordingly to where the file is
file_path2= "C:\\Users\\miyoh\\Git_Hub_Work\\Python_KLK6_data_analysis\\20240321_400CSF_selection_for_Toronto_reduced_FC_Corrected.xlsx" #clinical information of patients from collaborator
data= pd.read_excel(file_path, "Sheet2")
clinical_info = pd.read_excel(file_path2, "Sheet 1")
# %% [markdown]
# ## Preliminary checks
# %% [markdown]
# ### Prepare raw df (not clinical information dataframe)
# %%
data #check that it was read fine
#Result: it has 394 rows and 4 columns, but we should have 393 samples.
# %%
#the columns from df need to be renamed
data.columns = data.iloc[0] #the new column names are taken from row 1
data = data.drop(data.index[0]) #drop the first row, since we have now saved is as column names
# %%
data #check that it worked, new column names should be Value, ID, Group, Original_ID_ACE
#Results: 393 rows and 4 columns, Column headers are correct
# %%
#get info on data types and missing values
data.info()
#Results: no NAs, all types are object types
# %%
#make Value column a numeric value for proper subsequent analyses
# with this code, the other three columns remain intact, without being removed
data = data[["ID", "Group", "Value", "Original_ID_ACE"]].apply(pd.to_numeric, errors = "ignore")
# %%
#Define order of groups for later graphs
group_order = ['SCD', 'MCI non Progressor', 'MCI Progressor', 'AD']
# %% [markdown]
# ### Prepare clinical information
# %%
clinical_info
#correctly has 393 samples, 46 columns, but not all are needed
# %%
#Access column names, so that a selection can be made
clinical_info.columns
# %%
# Create another dataframe with picked columns
clinical_info_chosen = clinical_info[["Aliq", "MMSE_PL", "APOE", "Sex",
"F.Nacimiento", "FechaLCR", "Abeta_42_LCR",
"P_tau_LCR", "T_tau_LCR"]]
# %%
#check clinical_info_chosen
#this will show the top 5 and bottom 5 rows
clinical_info_chosen
#393 rows x9 columns
# %%
#get info on data types and missing values
clinical_info_chosen.info()
#Results: There is some missing data that will be dealt with later on
# %%
#turn F.Nacimiento (date of birth) into datetime()
clinical_info_chosen["F.Nacimiento"] = pd.to_datetime(clinical_info_chosen["F.Nacimiento"])
# %%
#Create a new column for Age = FechaLCR- F.Nacimiento
clinical_info_chosen['Age'] = [relativedelta(i,j).years for i,j in zip(clinical_info_chosen["FechaLCR"], clinical_info_chosen["F.Nacimiento"])]
# %%
#Rename Masculino = Male, Femenino = Female
clinical_info_chosen["Sex"] = clinical_info_chosen["Sex"].replace("Masculino", "Male")
clinical_info_chosen["Sex"] = clinical_info_chosen["Sex"].replace("Femenino", "Female")
# %%
#check results
print("The clinical data look like this:")
clinical_info_chosen
#Results: in column Sex, we should see Male, Female; in column Age, we should see the result of FechaLCR - F.Nacimiento as integer numbers.
# %% [markdown]
# ### Combine df with clinical_info_chosen
# %%
#Merge the dataframes with different column names
df = pd.merge(data, clinical_info_chosen, left_on= "Original_ID_ACE", right_on= "Aliq", how= "left")
print("\nMerged DataFrame (Left Join with different column names):")
print(df) #check result
#Result: 393 rows x 14 columns
# %% [markdown]
# ### Table 1
# %%
#Count percentage of female and male per Group
# Count and percentage by group
gender_percentage = df.groupby('Group')['Sex'].value_counts(normalize=True).mul(100).round(2) #round to 2 decimals
print("Percentage of each sex by group:")
print(gender_percentage)
# %%
#Count median of MMSE score per group
mmse_counts = df.groupby('Group')['MMSE_PL'].median().round(2) #round to 2 decimals
print("Median of the MMSE Score by group:")
print(mmse_counts)
# %%
#Count median of CSF Abeta 42 per group
Abeta_counts = df.groupby('Group')['Abeta_42_LCR'].median().round(2) #round to 2 decimals
print("Median of the Amyloid beta 1-42 by group:")
print(Abeta_counts)
# %%
#Count median of CSF p-tau per group
ptau_counts = df.groupby('Group')['P_tau_LCR'].median().round(2) #round to 2 decimals
print("Median of the p-tau by group:")
print(ptau_counts)
# %%
#Count median of CSF t-tau per group
ttau_counts = df.groupby('Group')['T_tau_LCR'].median().round(2) #round to 2 decimals
print("Median of the t-tau by group:")
print(ttau_counts)
# %%
#Count median of CSF hK6
hk6_counts = df.groupby("Group")['Value'].median().round(2)
print("hK6 median in each group:")
print(hk6_counts)
# %% [markdown]
# ### Modifications to df
# %%
# Change APOE status to correct versions with Greek lettering as per allele naming
df.loc[df['APOE'] == 'e2e3', 'APOE'] = 'ε2ε3'
df.loc[df['APOE'] == 'e2e4', 'APOE'] = 'ε2ε4'
df.loc[df['APOE'] == 'e3e3', 'APOE'] = 'ε3ε3'
df.loc[df['APOE'] == 'e3e4', 'APOE'] = 'ε3ε4'
df.loc[df['APOE'] == 'e4e4', 'APOE'] = 'ε4ε4'
# %% [markdown]
# ### Kruskal Wallis for all
# %%
#Perform Kruskal Wallis and Holm's correction
variables = ['Age', 'MMSE_PL', 'Abeta_42_LCR', 'P_tau_LCR', 'T_tau_LCR', 'Value'] # Variables to perform Kruskal Wallis on
# Create a results list to store all statistics
results = []
p_values = [] # Store p-values for later correction
print("Kruskal-Wallis Test Results for Multiple Variables")
print("=" * 60)
for var in variables:
print(f"\nVariable: {var}")
print("-" * 30)
# Extract data for each group
group_data = []
for group in group_order:
group_values = df[df['Group'] == group][var].dropna()
group_data.append(group_values)
print(f"Group {group}: n={len(group_values)}, median={group_values.median():.2f}")
# Perform Kruskal-Wallis test
h_stat, p_value = stats.kruskal(*group_data)
print(f"Kruskal-Wallis H-statistic: {h_stat:.3f}, p-value: {p_value:.4f}")
# Store results
result_dict = {
'Variable': var,
'H_statistic': h_stat,
'p_value': p_value,
'significant': p_value < 0.05
}
# Add group medians and sample sizes
for i, group in enumerate(group_order):
result_dict[f'Group_{group}_median'] = group_data[i].median()
result_dict[f'Group_{group}_n'] = len(group_data[i])
results.append(result_dict)
p_values.append(p_value) # Collect p-value for correction
# Apply Holm correction to all p-values AFTER collecting them all
rejected, holm_p_values, _, _ = multipletests(p_values, alpha=0.05, method='holm')
# Update results with Holm-corrected p-values
for i, result in enumerate(results):
result['Adjusted_p_value'] = holm_p_values[i]
result['Significant_after_correction'] = holm_p_values[i] < 0.05
# Convert to DataFrame for easier viewing
results_df = pd.DataFrame(results)
# Print summary of results
print("\n\nSummary with Holm-Bonferroni Correction:")
print("=" * 60)
for i, row in results_df.iterrows():
sig_symbol = "*" if row['Significant_after_correction'] else ""
print(f"{row['Variable']:15} p={row['p_value']:.4f} -> adj_p={row['Adjusted_p_value']:.4f} {sig_symbol}")
# %% [markdown]
# ## Check normality, variance (homoscedasticity) of hK6 data
# %% [markdown]
# To see if the data are normally distributed, it is possible to run two statistical tests:
# - Shapiro-Wilk Test, which is suitable for smaller sample sizes (typically <50), where the null hypothesis is that the data are normally distributed.
# - D'Agostino's K-squared test (normaltest), which is based on skewness and kurtosis and is generally suitable for larger sample sizes. Similar to Shapiro-Wilk, a low p-value suggests non-normalilty.
# %%
#Perform D'Agostino K^2 test
groups = df['Group'].unique() #group by Group column
for group_name in groups:
group_data = df[df["Group"] == group_name]["Value"] #group the values based on the Groups
statistic, p_value = stats.normaltest(group_data) #if we wanted to do Shapiro-Wilk, we would run a shapiro()
print(f"Normality test for Group {group_name}:")
print(f" D'Agostino K^2 Statistic: {statistic:.4f}")
print(f" P-value: {p_value:.4f}")
alpha = 0.05 #significance level
if p_value > alpha:
print(f"Result: Data for Group {group_name} appears to be normally distributed (fail to reject H0)")
else:
print(f"Result: Data for group {group_name} does not appear to be normally distributed (reject H0)")
print("-" * 30)
# %%
#Alternatively, check normal distribution using QQ plot
#It helps to visualize the distribution
fix, axes= plt.subplots(2,2, figsize=(8,6), dpi=600)
axes=axes.flatten() #Flatten the 2x2 array of axes for easy indexing
for i, group_name in enumerate(groups):
if i <len(axes): #this ensures that we don't exceed the number of subplots
group_data = df[df['Group'] == group_name]['Value'].dropna()
stats.probplot(group_data, dist= 'norm', plot=axes[i])
axes[i].set_title(f"Q-Q Plot for {group_name} Group \n (n= {len(group_data)})")
plt.tight_layout()
plt.suptitle("Q-Q Plots for Normality Assessment", y= 1.02, fontsize= 14)
plt.show()
# %%
#Since our data is not normally distributed, we log2 transform our values
#in a new column, which will be used from now on for any parametric tests chosen
df["log2_Value"] = np.log2(df["Value"])
# %%
#Check for normality again using the log2 values
for group_name in groups:
group_data = df[df["Group"] == group_name]["log2_Value"]
statistic, p_value = stats.normaltest(group_data)
print(f"Normality test for Group {group_name}:")
print(f" D'Agostino Statistic: {statistic:.4f}")
print(f" P-value: {p_value:.4f}")
alpha = 0.05 #significance level
if p_value > alpha:
print(f"Result: Data for Group {group_name} appears to be normally distributed (fail to reject H0)")
else:
print(f"Result: Data for group {group_name} does not appear to be normally distributed (reject H0)")
print("-" * 30)
#Result: Our hK6 data is now normally distributed
# %%
# I will also check normality for age, because we will have to adjust for it
for group_name in groups:
group_data = df[df['Group'] == group_name]['Age']
statistic, p_value = stats.normaltest(group_data)
print(f'Normality test for group {group_name}:')
print(f"D'Agostino Statistic: {statistic:.4f}")
print(f"P-value: {p_value}")
alpha=0.05
if p_value > alpha:
print(f"Result: Data for Group {group_name} appears to be normally distributed (fail to reject H0)")
else:
print(f"Result: Data for group {group_name} does not appear to be normally distributed (reject H0)")
print("-" * 30)
# %% [markdown]
# Note: ANOVA is more powerful than Kruskal-Wallis, since Kruskal-Wallis compares the ranks, and some information may be lost. When we have ordinal data (such as MMSE Score)
# we have to use Kruskal-Wallis.
# %% [markdown]
# Prior to ANOVA, we have to check for variances, since they have to be equal (we have to assume that the variability within each group being compared is roughly the same).
# %% [markdown]
# To do this, we can do the Bartlett's test, or Levene's test.
# We have used Levene's test which checks for homogeneity of variances but is less sensitive to departures from normality than Bartlett's test, making it more robust.
# %%
#Create 4 groups for stat analysis based on Group column
group_nonprog = df[df["Group"] == "MCI non Progressor"][['Group','log2_Value',"Value", "MMSE_PL", "APOE", "Sex","Age","Abeta_42_LCR", "P_tau_LCR", "T_tau_LCR"]]
group_prog = df[df["Group"] == "MCI Progressor"][['Group','log2_Value',"Value", "MMSE_PL", "APOE", "Sex","Age","Abeta_42_LCR", "P_tau_LCR", "T_tau_LCR"]]
group_SCD= df[df["Group"] == "SCD"][['Group', 'log2_Value',"Value", "MMSE_PL", "APOE", "Sex","Age","Abeta_42_LCR", "P_tau_LCR", "T_tau_LCR"]]
group_AD = df[df["Group"] == "AD"][['Group', 'log2_Value',"Value", "MMSE_PL", "APOE", "Sex","Age","Abeta_42_LCR", "P_tau_LCR", "T_tau_LCR"]]
# %%
#Check for NAs for the Levene's test
print(f"nonprog: {group_nonprog['log2_Value'].isna().sum()}")
print(f"prog: {group_prog['log2_Value'].isna().sum()}")
print(f"SCD: {group_SCD['log2_Value'].isna().sum()}")
print(f"AD: {group_AD['log2_Value'].isna().sum()}")
# %%
# Check variances
print("Group variances:")
print(f"nonprog: {group_nonprog['log2_Value'].var()}")
print(f"prog: {group_prog['log2_Value'].var()}")
print(f"SCD: {group_SCD['log2_Value'].var()}")
print(f"AD: {group_AD['log2_Value'].var()}")
# Check if any group has zero variance (all values identical)
print("All values same in groups?")
print(f"nonprog: {group_nonprog['log2_Value'].nunique() == 1}")
print(f"prog: {group_prog['log2_Value'].nunique() == 1}")
print(f"SCD: {group_SCD['log2_Value'].nunique() == 1}")
print(f"AD: {group_AD['log2_Value'].nunique() == 1}")
# %%
#perform Levene's test
l_statistic, p_value_l = stats.levene(group_nonprog['log2_Value'], group_prog['log2_Value'], group_SCD['log2_Value'], group_AD['log2_Value'])
print(f"Levene's Test Statistic: {l_statistic}")
print(f"Levene's Test p-value: {p_value_l}")
# Interpretation: If p_value < significance level (e.g., 0.05), reject the null hypothesis
# and conclude that variances are not equal.
print("-"*30) #divider
alpha = 0.05 #significance level
if p_value_l > alpha:
print("The variances are equal (fail to reject the H0)")
else:
print("The variances are not equal (reject the H0)")
#Result: The variances are equal
# %% [markdown]
# ## Performing ANCOVA for the hK6 data
# %% [markdown]
# We have to perform the statistics, in this case a rank-based ANCOVA test (pnon-arametric), using the log2 transformed data
# %%
# 2. Check distribution before/after transformation
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
sns.histplot(df['Value'], ax=ax1, kde=True)
ax1.set_title('Raw hK6 Distribution')
sns.histplot(df['log2_Value'], ax=ax2, kde=True)
ax2.set_title('Log2-Transformed hK6 Distribution')
plt.show()
# 3. Perform ANCOVA on transformed data
model = smf.ols('log2_Value ~ Age + C(Group)', data=df).fit()
ancova_table = sm.stats.anova_lm(model, typ=2)
print("ANCOVA Results:")
print(ancova_table)
# 4. Check model assumptions
residuals = model.resid
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
stats.probplot(residuals, dist="norm", plot=axes[0])
axes[0].set_title('Q-Q Plot of Residuals')
axes[1].scatter(model.fittedvalues, residuals)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_xlabel('Fitted values')
axes[1].set_ylabel('Residuals')
axes[1].set_title('Residuals vs Fitted')
plt.show()
df['hK6_rank'] = df['log2_Value'].rank()
# ANCOVA model: Ranked_KLK6 ~ Age + Diagnosis
model_rank = smf.ols('hK6_rank ~ Age + C(Group)', data=df).fit()
print("ANCOVA Results:")
print(sm.stats.anova_lm(model_rank, typ=2))
# %% [markdown]
# Based on the above results, we should stick to the parametric ANCOVA.
# %% [markdown]
# In order to see which groups are significantly different from each other, a post-hoc test is used. Tukey's is a common choice, which we have used here.
# %%
tukey_result = pairwise_tukeyhsd(endog=df['log2_Value'], groups = df['Group'], alpha = 0.05)
print(tukey_result)
p_values_tukey = tukey_result.pvalues
# %%
p_values_tukey #check to see what the different p_values are to use later
# %% [markdown]
# ### Function for hK6 plot asterisks
# %%
# Function to add significance bars with custom symbols, this will be used in the creation of the graphs for hK6 later on.
'''
Parameters
ax= name of plot that has already been defined
group1, group2= the two groups that we are comparing (we don't need to mention the dataframe, since it has been typed in ax)
p_value= the p_value from the Tukey test, include location in [], e.g. p_values_tukey[0] for AD vs MCI_conver_toDem
x1, x2= location of boxplot for each group compared, e.g. MCI_NonProgressor is in x1= 0, while MCI_conver_toDem is in x1= 1.
y= height of line based on the graph.
'''
def add_sig_bars(ax, group1, group2, p_value, x1, x2, y):
"""Add significance bars with custom symbols and sizes based on p-value."""
y_offset = 0.2 # Offset for the bar
ax.plot([x1, x1, x2, x2], [y, y + y_offset, y + y_offset, y], color='black')
# Annotate with different symbols based on p-value
if p_value < 0.001:
ax.text((x1 + x2) / 2, y + y_offset + 0.1, '***', fontsize=12, ha='center', color='red')
elif p_value < 0.01:
ax.text((x1 + x2) / 2, y + y_offset + 0.1, '**', fontsize=12, ha='center', color='orange')
elif p_value < 0.05:
ax.text((x1 + x2) / 2, y + y_offset + 0.1, '*', fontsize=12, ha='center', color='blue')
# %% [markdown]
# ## hK6 plots
# NOTE: MCI_conver_toDem corresponds to MCI Progressor in the later analysis
# %% [markdown]
# ### Create plots with raw data
# %%
#Scatterplot with boxplot
plt.figure(figsize=(10,6))
# Set the style to include a grid
sns.set_style("whitegrid")
ax1 = sns.boxplot(y = "Value", x= 'Group',palette= "RdBu", data= df,order= group_order)
sns.stripplot(y = "Value", x="Group", ax= ax1, data= df, color = "black")
plt.title("Comparison of hK6 in the AD continuum", fontsize = 17)
# %%
#Boxplot without the scatterplot
plt.figure(figsize=(10,6))
ax = sns.boxplot(y = df["Value"], x= df["Group"], palette= "RdBu", order= group_order)
# Set the style to include a grid
sns.set_style("whitegrid")
plt.title("Comparison of hK6 in the AD continuum", fontsize = 17)
# %% [markdown]
# ### Create plots with log2 transformed data
# %%
#Scatterplot with boxplot
plt.figure(figsize=(10,6))
# Set the style to include a grid
sns.set_style("whitegrid")
ax1 = sns.boxplot(y = df["log2_Value"], x= df["Group"], hue= "Group",palette= "RdBu", data= df,
hue_order= group_order, order = group_order)
sns.stripplot(y = df["log2_Value"], x= df["Group"], ax= ax1, data= df, color = ".3")
#manually check the sig p values
add_sig_bars(ax1, "MCI non Progessor", "MCI Progressor",p_values_tukey[3],0 ,1, 9.5)
add_sig_bars(ax1, "AD", "MCI non Progessor",p_values_tukey[1],0 ,3, 10)
add_sig_bars(ax1, "AD", "SCD",p_values_tukey[2],2 ,3, 10.4)
add_sig_bars(ax1, "SCD", "MCI Progressor",p_values_tukey[4],1 ,2, 9.4)
plt.title("Comparison of hK6 in the AD continuum", fontsize = 17)
plt.ylim(top=11.0)
# %%
#Boxplot without scatterplot
plt.figure(figsize=(10,6))
ax = sns.boxplot(y = df["log2_Value"], x= df["Group"], palette= "RdBu",
hue_order= group_order, order = group_order)
# Set the style to include a grid
sns.set_style("whitegrid")
#manually check the sig p values
add_sig_bars(ax, "MCI non Progessor", "MCI Progressor",p_values_tukey[3],0 ,1, 9.5)
add_sig_bars(ax, "AD", "MCI non Progessor",p_values_tukey[1],0 ,3, 10)
add_sig_bars(ax, "AD", "SCD",p_values_tukey[2],2 ,3, 10.4)
add_sig_bars(ax, "SCD", "MCI Progressor",p_values_tukey[4],1 ,2, 9.4)
plt.title("Comparison of hK6 in the AD continuum", fontsize = 17)
plt.ylim(top=11.0)
# %% [markdown]
# ## Function for correlation analysis
# %%
#CREATE A FUNCTION FOR SPEARMAN CORRELATION FIGURES AND RANK-TRANSFORMED FIGURES
'''
Parameters
x,y= variable names or arrays
data= DataFrame that contains the variables (optional)
title= Plot title (Optional)
xlabel, ylabel= Axis labels (optional)
group= Grouping variable for color coding (optional)
group_order= Order of groups for coloring (optional)
'''
def plot_spearman_correlation(x, y, data=None, title=None, xlabel=None, ylabel=None, group=None, group_order= None):
if data is not None:
x_data = data[x] #the x value from the dataframe
y_data = data[y] #the y data from the dataframe
x_col_name= x #the x column name
y_col_name= y #the y column name
if group is not None:
group_data = data[group] #if group data is provided, then use it
else:
x_data = x
y_data = y
x_col_name= "X Variable"
y_col_name= "Y Variable"
group_data= group
# use provided labels or use the default ones as above
xlabel = xlabel if xlabel is not None else x_col_name
ylabel = ylabel if ylabel is not None else y_col_name
# Calculate spearman correlation
corr, p_value = stats.spearmanr(x_data, y_data)
# Create figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6), dpi=600)
#color coding based on the groups
if group is not None:
#Get unique groups and then colors
if group_order is not None:
unique_groups= group_order
else:
unique_groups = group_data.unique()
n_groups= len(unique_groups) #number of groups
colors= sns.color_palette("RdBu", n_groups) #choose the colour for each of the groups
#Plot each group with a different color
for i, grp in enumerate(unique_groups): #grp=group
mask = group_data == grp
ax1.scatter(x_data[mask], y_data[mask], alpha= 0.6, s=20,
color= colors[i], label= str(grp))
ax2.scatter(x_data.rank()[mask], y_data.rank()[mask], alpha=0.6,
s=20, color= colors[i], label= str(grp))
#add legend code here if needed
else:
#code without groups
ax1.scatter(x_data, y_data, alpha=0.6, s=20)
ax2.scatter(x_data.rank(), y_data.rank(), alpha=0.6, s=20)
# Add regression lines to both plots
# For original data
slope_orig, intercept_orig, r_value_orig, p_value_orig, std_err_orig = stats.linregress(x_data, y_data)
x_range_orig = np.linspace(x_data.min(), x_data.max(), 100)
y_pred_orig = slope_orig * x_range_orig + intercept_orig
ax1.plot(x_range_orig, y_pred_orig, '-',color='black', linewidth=2, alpha=0.8,
label=f'Slope: {slope_orig:.2f}')
# For rank-transformed data
x_rank = x_data.rank()
y_rank = y_data.rank()
slope_rank, intercept_rank, r_value_rank, p_value_rank, std_err_rank = stats.linregress(x_rank, y_rank)
x_range_rank = np.linspace(x_rank.min(), x_rank.max(), 100)
y_pred_rank = slope_rank * x_range_rank + intercept_rank
ax2.plot(x_range_rank, y_pred_rank, '-',color='black', linewidth=2, alpha=0.8,
label='Slope')
# Original data scatter plot
ax1.set_xlabel(xlabel)
ax1.set_ylabel(ylabel)
ax1.set_title('Original Data')
ax1.grid(True, alpha=0.3)
# Rank-transformed data
x_rank = x_data.rank()
y_rank = y_data.rank()
ax2.set_xlabel(f'{xlabel} Rank')
ax2.set_ylabel(f'{ylabel} Rank')
ax2.set_title('Rank-Transformed Data')
ax2.grid(True, alpha=0.3)
# Add perfect correlation line to rank plot (1:1 line) depending on whether it is positive or negative
if slope_rank >0:
min_rank = min(x_rank.min(), y_rank.min())
max_rank = max(x_rank.max(), y_rank.max())
ax2.plot([min_rank, max_rank], [min_rank, max_rank], '--',color= 'purple', alpha=0.7,
label= 'Perfect pos slope: 1')
ax2.legend(bbox_to_anchor=(-0.1, 1)) #add legend
else:
min_rank = min(x_rank.min(), y_rank.min())
max_rank = max(x_rank.max(), y_rank.max())
ax2.plot([min_rank, max_rank], [max_rank, min_rank], '--',color= 'purple', alpha=0.7,
label= 'Perfect neg slope: -1')
ax2.legend(bbox_to_anchor=(-0.1, 1)) #add legend
# Add correlation info to both plots
for ax in [ax1, ax2]:
ax.text(0.05, 0.95, f'ρ = {corr:.2f}\nraw p = {p_value:.3f}',
transform=ax.transAxes,
bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.9))
plt.suptitle(f'{title} (Spearman ρ = {corr:.2f}, p= {p_value:.3f})', fontsize=14)
plt.tight_layout()
plt.savefig(f'{title}.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ## Correlation of hK6 with Age, Sex, MMSE score, APOE status, Abeta, P tau and T tau (combined data)
# %% [markdown]
# - hK6 protein concentration is a continuous variable, Age is a continous variable, Sex is categorical (Binary), MMSE score is continous, and APOE status is categorical (Ordinal/Nominal).
# - When comparing two continous variables, we performed a spearman correlation.
# - When comparing continous to categorical (Binary) we performed the Point-Biserial correlation.
# - When comparing continous to categorical (Ordinal/Nominal) we performed an Kruskal-Wallis.
# %%
df #check how it looks
#Result: 393 rows x 15 columns
# %%
#Check descriptive stats
print(df.describe(include='all'))
# %%
#Check specificically for categorical variables
print("\nValue counts for categorical variables:")
print("Sex:\n", df['Sex'].value_counts())
print("\nAPOE_Status:\n", df['APOE'].value_counts())
# %%
### Check correlations of all variables with age
corr_matrix = df[['Age', 'Value', 'MMSE_PL', 'Abeta_42_LCR', 'P_tau_LCR', 'T_tau_LCR']].corr(method = 'spearman')
print("correlation matrix with Age: ")
print(corr_matrix)
# %% [markdown]
# There is a positive correlation with p-tau and t-tau and hK6 concentration, and a negative correlation with MMSE score and Abeta 42.
# %% [markdown]
# Thus, we should adjust for age through a partial correlation.
# %%
# Create a binary variable for sex (0/1)
df['Sex_numeric'] = df['Sex'].map({'Male': 0, 'Female': 1}) # or vice versa
# %% [markdown]
# ### Correlation with Age
# %%
#Correlation with age
print("=" *50)
print("CORRELATION WITH AGE")
print("="*50)
# Check normality using QQ plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
stats.probplot(df['Value'].dropna(), dist="norm", plot=ax1)
ax1.set_title('QQ Plot - hK6 Levels') #all diagnoses combined
stats.probplot(df['Age'].dropna(), dist="norm", plot=ax2)
ax2.set_title('QQ Plot - Age')
plt.tight_layout()
plt.show()
# %%
##Correlation with Age
print('='*50)
print('CORRELATION WITH AGE')
print('='*50)
#Fix dropna issues (the x and y variables are not the same length)
age_temp_df = df[["log2_Value","Value", 'Age', 'Sex_numeric', 'Group']].dropna()
#Spearman correlation because Abeta 42 is a continous variable
spearman_corr_age, spearman_p_age = stats.spearmanr(age_temp_df['Value'], age_temp_df['Age'])
print(f"Spearman Correlation: ρ = {spearman_corr_age:.3f}, p = {spearman_p_age:.4f}")
print(f'The sample size is: {len(age_temp_df)}')
#Adjust for sex
print('='*50)
print("CORRELATION WITH AGE ( ADJUSTED FOR SEX)")
print("="*50)
partial_spearman_age = pg.partial_corr(data= age_temp_df, x= 'Age', y='Value', covar= 'Sex_numeric', method= 'spearman')
print('Partial correlation (controlling for sex):')
print(f"ρ = {partial_spearman_age['r'].iloc[0]:.3f}, p = {partial_spearman_age['p-val'].iloc[0]:.4f}")
print('*'*50)
# Compare unadjusted vs adjusted results
print(f"\nComparison:")
print(f"Unadjusted p-value (Mann-Whitney): {spearman_p_age:.4f}")
print(f"Age-adjusted p-value (ANCOVA): {partial_spearman_age['p-val'].iloc[0]:.4f}")
if abs(spearman_p_age - partial_spearman_age['p-val'].iloc[0]) > 0.05:
print("Note: Adjustment for age meaningfully changed the results")
else:
print("Note: Results were robust to age adjustment")
#Create a scatterplot
plt.figure(figsize=(8,6), dpi= 600)
sns.scatterplot(data= age_temp_df, x= "Age", y= "Value",
hue= 'Group', palette= 'RdBu', hue_order= group_order)
plt.title(f"hK6 vs. Age (Spearman ρ = {spearman_corr_age:.3f}, p = {spearman_p_age:.4f})")
plt.xlabel("Age (years)")
plt.ylabel("hK6 concentration (ng/mL)")
plt.savefig('comb_hK6vsAge.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Spearman correlation when combining all data
plot_spearman_correlation("Age", "Value", data= df, title= "hK6 vs. Age (Spearman Correlation)",
xlabel= "Age (years)",ylabel= "hK6 concentration (ng/mL)", group= 'Group', group_order= group_order)
# %% [markdown]
# ### Correlation with Sex
# %%
#Correlation with Sex
print('='*50)
print('CORRELATION WITH SEX')
print('='*50)
# Check normality within each group using Shapiro-Wilk
for sex_group in df['Sex'].unique():
group_data = df[df['Sex'] == sex_group]['Value'].dropna()
stat, p = stats.shapiro(group_data)
print(f"Shapiro-Wilk for {sex_group}: p = {p:.4f}")
#Use Mann-Whitney test, non-normal data
male_data = df[df['Sex']== 'Male']['Value'].dropna()
female_data = df[df['Sex']== 'Female']['Value'].dropna()
stat_sex, p_value_sex = stats.mannwhitneyu(male_data, female_data, alternative= 'two-sided')
print(f"\nMann-Whitney U test: U = {stat_sex}, p = {p_value_sex:.4f}")
print('*'*50)
print('='*50)
print('CORRELATION WITH SEX (ADJUSTED FOR AGE)')
print('='*50)
print("\n2. Partial Correlation Approach:")
partial_corr_sex = pg.partial_corr(data=df, x='Sex_numeric', y='Value', covar='Age', method='spearman')
print(f"Partial correlation (Sex vs Value, controlling for Age):")
print(f"ρ = {partial_corr_sex['r'].iloc[0]:.3f}, p = {partial_corr_sex['p-val'].iloc[0]:.4f}")
# Compare unadjusted vs adjusted results
print(f"\nComparison:")
print(f"Unadjusted p-value (Mann-Whitney): {p_value_sex:.4f}")
print(f"Age-adjusted p-value (ANCOVA): {partial_corr_sex['p-val'].iloc[0]:.4f}")
if abs(p_value_sex - partial_corr_sex['p-val'].iloc[0]) > 0.05:
print("Note: Adjustment for age and meaningfully changed the results")
else:
print("Note: Results were robust to age adjustment")
# %%
# Create a box plot
plt.figure(figsize=(8, 6), dpi=600)
sns.boxplot(data=df, x='Sex', y='Value',boxprops= dict(alpha=0.6), hue= "Sex", palette= "RdBu")
sns.stripplot(data=df, x='Sex', y='Value', color='black', alpha=0.5, size=4)
plt.xlabel("Sex")
plt.ylabel("hK6 concentration (ng/mL)")
plt.title(f'hK6 Levels by Sex (Mann-Whitney p = {p_value_sex:.3f})')
plt.savefig('comb_hK6vsSex.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correlation with MMSE score
# %%
##Correlation with MMSE Score
print('='*50)
print('CORRELATION WITH MMSE Score')
print('='*50)
#Fix dropna issues (the x and y variables are not the same length)
mmse_temp_df = df[["log2_Value","Value", 'MMSE_PL', 'Sex_numeric', 'Age', 'Group']].dropna()
#Spearman correlation because Abeta 42 is a continous variable
spearman_corr_mmse, spearman_p_mmse = stats.spearmanr(mmse_temp_df['Value'], mmse_temp_df['MMSE_PL'])
print(f"Spearman Correlation: ρ = {spearman_corr_mmse:.3f}, p = {spearman_p_mmse:.4f}")
print(f'The sample size is: {len(mmse_temp_df)}')
#Adjust for sex
print('='*50)
print("CORRELATION WITH MMSE Score ( ADJUSTED FOR SEX AND AGE)")
print("="*50)
partial_spearman_mmse = pg.partial_corr(data= mmse_temp_df, x= 'MMSE_PL', y='Value',
covar= ['Sex_numeric', 'Age'], method= 'spearman')
print('Partial correlation (controlling for sex and age):')
print(f"ρ = {partial_spearman_mmse['r'].iloc[0]:.3f}, p = {partial_spearman_mmse['p-val'].iloc[0]:.4f}")
print('*'*50)
# Compare unadjusted vs adjusted results
print(f"\nComparison:")
print(f"Unadjusted p-value (Mann-Whitney): {spearman_p_mmse:.4f}")
print(f"Age-adjusted p-value (ANCOVA): {partial_spearman_mmse['p-val'].iloc[0]:.4f}")
if abs(spearman_p_mmse - partial_spearman_mmse['p-val'].iloc[0]) > 0.05:
print("Note: Adjustment for age and sex meaningfully changed the results")
else:
print("Note: Results were robust to age and sex adjustment")
#Create a scatterplot
plt.figure(figsize=(8,6), dpi= 600)
sns.scatterplot(data= mmse_temp_df, x= "MMSE_PL", y= "Value",
hue= 'Group', palette= 'RdBu', hue_order= group_order)
plt.title(f"hK6 vs. MMSE Score (Spearman ρ = {spearman_corr_mmse:.3f}, p = {spearman_p_mmse:.4f})")
plt.xlabel("MMSE Score (1-30 scale)")
plt.ylabel("hK6 concentration (ng/mL)")
plt.savefig('comb_hK6vsMMSE.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
## Repeat analysis by checking the covariates separately
# Correlation with MMSE Score
print('='*50)
print('CORRELATION WITH MMSE Score')
print('='*50)
# Fix dropna issues
mmse_temp_df = df[["log2_Value","Value", 'MMSE_PL', 'Sex_numeric', 'Age', 'Group']].dropna()
# Spearman correlation
spearman_corr_mmse, spearman_p_mmse = stats.spearmanr(mmse_temp_df['Value'], mmse_temp_df['MMSE_PL'])
print(f"Spearman Correlation: ρ = {spearman_corr_mmse:.3f}, p = {spearman_p_mmse:.4f}")
print(f'The sample size is: {len(mmse_temp_df)}')
# ===========================================================================
# STEP 1: TEST EACH COVARIATE SEPARATELY
# ===========================================================================
print('='*50)
print("IDENTIFYING WHICH COVARIATE CAUSES SIGNIFICANCE LOSS")
print("="*50)
# 1. Adjust for SEX only
partial_sex_only = pg.partial_corr(data=mmse_temp_df, x='MMSE_PL', y='Value',
covar=['Sex_numeric'], method='spearman')
# 2. Adjust for AGE only
partial_age_only = pg.partial_corr(data=mmse_temp_df, x='MMSE_PL', y='Value',
covar=['Age'], method='spearman')
# 3. Adjust for BOTH (your original)
partial_both = pg.partial_corr(data=mmse_temp_df, x='MMSE_PL', y='Value',
covar=['Sex_numeric', 'Age'], method='spearman')
print(f"\nUnadjusted: ρ = {spearman_corr_mmse:.3f}, p = {spearman_p_mmse:.4f}")
print(f"Adjusted for SEX only: ρ = {partial_sex_only['r'].iloc[0]:.3f}, p = {partial_sex_only['p-val'].iloc[0]:.4f}")
print(f"Adjusted for AGE only: ρ = {partial_age_only['r'].iloc[0]:.3f}, p = {partial_age_only['p-val'].iloc[0]:.4f}")
print(f"Adjusted for BOTH: ρ = {partial_both['r'].iloc[0]:.3f}, p = {partial_both['p-val'].iloc[0]:.4f}")
# ===========================================================================
# STEP 2: DETERMINE WHICH COVARIATE IS RESPONSIBLE
# ===========================================================================
print('\n' + '*'*50)
print("ANALYSIS OF SIGNIFICANCE CHANGE")
print('*'*50)
# Check if significance was lost
if spearman_p_mmse < 0.05 and partial_both['p-val'].iloc[0] >= 0.05:
print("✓ Significance was LOST after adjustment")
# Check which adjustment caused the change
sex_change = abs(spearman_p_mmse - partial_sex_only['p-val'].iloc[0])
age_change = abs(spearman_p_mmse - partial_age_only['p-val'].iloc[0])
print(f"Change due to SEX adjustment: Δp = {sex_change:.4f}")
print(f"Change due to AGE adjustment: Δp = {age_change:.4f}")
if sex_change > age_change and partial_sex_only['p-val'].iloc[0] >= 0.05:
print("→ PRIMARY CULPRIT: SEX (adjusting for sex alone removes significance)")
elif age_change > sex_change and partial_age_only['p-val'].iloc[0] >= 0.05:
print("→ PRIMARY CULPRIT: AGE (adjusting for age alone removes significance)")
elif partial_sex_only['p-val'].iloc[0] < 0.05 and partial_age_only['p-val'].iloc[0] >= 0.05:
print("→ PRIMARY CULPRIT: AGE (sex adjustment preserves significance)")
elif partial_age_only['p-val'].iloc[0] < 0.05 and partial_sex_only['p-val'].iloc[0] >= 0.05:
print("→ PRIMARY CULPRIT: SEX (age adjustment preserves significance)")
else:
print("→ COMBINED EFFECT: Both age and sex contribute to significance loss")
elif spearman_p_mmse >= 0.05 and partial_both['p-val'].iloc[0] < 0.05:
print("✓ Significance was GAINED after adjustment")
# Similar logic for gained significance...
else:
print("✓ No major change in significance after adjustment")
# ===========================================================================
# STEP 3: CHECK FOR CONFOUNDING RELATIONSHIPS
# ===========================================================================
print('\n' + '='*50)
print("CONFOUNDING ANALYSIS")
print("="*50)
# Check correlations between covariates and main variables
print("Correlations to identify potential confounding:")
print(f"MMSE vs Age: ρ = {stats.spearmanr(mmse_temp_df['MMSE_PL'], mmse_temp_df['Age'])[0]:.3f}")
print(f"MMSE vs Sex: ρ = {stats.spearmanr(mmse_temp_df['MMSE_PL'], mmse_temp_df['Sex_numeric'])[0]:.3f}")
print(f"Value vs Age: ρ = {stats.spearmanr(mmse_temp_df['Value'], mmse_temp_df['Age'])[0]:.3f}")
print(f"Value vs Sex: ρ = {stats.spearmanr(mmse_temp_df['Value'], mmse_temp_df['Sex_numeric'])[0]:.3f}")
# ===========================================================================
# STEP 4: VISUALIZE THE RELATIONSHIPS
# ===========================================================================
# Create faceted plots to see the relationships
fig, axes = plt.subplots(2, 2, figsize=(12, 10), dpi=600)
# Plot 1: Original relationship
sns.scatterplot(data=mmse_temp_df, x="MMSE_PL", y="Value", hue='Group',
palette='RdBu', ax=axes[0,0])
axes[0,0].set_title(f"Original: ρ = {spearman_corr_mmse:.3f}, p = {spearman_p_mmse:.4f}")
# Plot 2: Stratified by sex
sns.scatterplot(data=mmse_temp_df, x="MMSE_PL", y="Value", hue='Sex_numeric',
palette='viridis', ax=axes[0,1])
axes[0,1].set_title("Colored by Sex")
# Plot 3: Colored by age (categorized)
mmse_temp_df['Age_group'] = pd.cut(mmse_temp_df['Age'], bins=3)
sns.scatterplot(data=mmse_temp_df, x="MMSE_PL", y="Value", hue='Age_group',
palette='plasma', ax=axes[1,0])
axes[1,0].set_title("Colored by Age Group")
# Plot 4: Residuals after adjusting for age and sex (optional)
from sklearn.linear_model import LinearRegression
X_adjust = mmse_temp_df[['Age', 'Sex_numeric']]
model = LinearRegression().fit(X_adjust, mmse_temp_df['Value'])
residuals = mmse_temp_df['Value'] - model.predict(X_adjust)
axes[1,1].scatter(mmse_temp_df['MMSE_PL'], residuals, c=mmse_temp_df['Age'], cmap='plasma')
axes[1,1].set_xlabel("MMSE Score")
axes[1,1].set_ylabel("Value (adjusted for Age/Sex)")
axes[1,1].set_title("Relationship with MMSE after adjusting for Age/Sex")
plt.tight_layout()
plt.savefig('MMSE_confounding_analysis.jpg', dpi=600, bbox_inches='tight')
plt.show()
print('\n' + '='*50)
print("FINAL CONCLUSION")
print("="*50)
print("Check the individual adjustments above to see which covariate")
print("caused the greatest change in p-value and effect size.")
# %%
#Correlation with MMSE score
plot_spearman_correlation("MMSE_PL", "Value", data= mmse_temp_df, xlabel= "MMSE Score", ylabel =
"hK6 concentration (ng/mL)",
title= "hK6 vs. MMSE Score (Spearman correlation)",
group= 'Group',
group_order=group_order)
# %% [markdown]
# ### Correlation with APOE status
# %%
#Correlation with APOE status, performing ANOVA on the log2 hK6 values
print('='*50)
print("CORRELATION WITH APOE STATUS")
print('='*50)
#remove NAs like in the MMSE analysis
apoe_temp_df = df[["log2_Value","Value", 'APOE','Age', 'Sex_numeric', 'Group']].dropna()
# Prepare data for ANOVA - create a list of arrays for each group
apoe_groups = []
for status in apoe_temp_df['APOE'].unique():
group_data = apoe_temp_df[apoe_temp_df['APOE'] == status]['log2_Value']
apoe_groups.append(group_data)
print(f"APOE {status}: n = {len(group_data)}, mean = {np.mean(group_data):.2f} ± {np.std(group_data):.2f}")
# Perform one-way ANOVA
f_stat_apoe, p_value_apoe = f_oneway(*apoe_groups)
print(f"\nOne-way ANOVA: F = {f_stat_apoe:.3f}, p = {p_value_apoe:.4f}")
# If ANOVA is significant (p < 0.05), perform Tukey's HSD post-hoc test
if p_value_apoe < 0.05:
print("\nPerforming Tukey's HSD post-hoc test:")
# Prepare data for Tukey test
tukey_data = apoe_temp_df[['log2_Value', 'APOE']].dropna()
tukey_results = pairwise_tukeyhsd(tukey_data['log2_Value'], tukey_data['APOE'])
print(tukey_results)
# Plot the results
tukey_results.plot_simultaneous()
plt.title("Tukey HSD Post-hoc Test - Confidence Intervals")
plt.show()
else:
print("ANOVA not significant - no post-hoc tests needed.")
# ADJUSTED ANALYSIS - ANCOVA approach
print('='*50)
print('APOE ANALYSIS (ADJUSTED FOR AGE AND SEX)')
print('='*50)
# First, describe the groups
apoe_groups = []
for status in sorted(apoe_temp_df['APOE'].unique()):
group_data = apoe_temp_df[apoe_temp_df['APOE'] == status]['log2_Value'].dropna()
apoe_groups.append(group_data)
print(f"APOE {status}: n = {len(group_data)}, mean = {np.mean(group_data):.2f} ± {np.std(group_data):.2f}")
# Perform ANCOVA (correct approach for categorical variables)
ancova_result = pg.ancova(data=apoe_temp_df, dv='log2_Value', between='APOE', covar=['Age', 'Sex_numeric'])
print('\nANCOVA Results (adjusted for Age and Sex):')
print(ancova_result)
# Post-hoc tests if ANCOVA is significant
if ancova_result['p-unc'].iloc[0] < 0.05:
print('\nPost-hoc pairwise comparisons:')
posthoc = pg.pairwise_ancova(data=apoe_temp_df, dv='log2_Value', between='APOE', covar=['Age', 'Sex_numeric'])
print(posthoc)
# %%
# Create a box plot, plot raw values instead of log2_values
plt.figure(figsize=(10, 6), dpi=600)
sns.boxplot(data=apoe_temp_df, x='APOE', y='Value', hue = "APOE",
palette= "viridis", boxprops = dict(alpha= 0.6))
sns.stripplot(data=apoe_temp_df, x='APOE', y='Value', color='black', alpha=0.5, size=4)
plt.title(f'hK6 Levels by APOE Status (ANOVA p = {p_value_apoe:.3f})')
plt.xticks(rotation=45)
plt.xlabel("APOE Status")
plt.ylabel("hK6 concentration (ng/mL)")
plt.tight_layout()
plt.savefig('comb_hK6vsAPOE.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correlationn with Abeta 42
# %%
#Correlation with Abeta 42
print('='*50)
print('CORRELATION WITH ABETA 42')
print('='*50)
#Fix dropna issues (the x and y variables are not the same length)
abeta_temp_df = df[["log2_Value","Value", 'Abeta_42_LCR','Sex_numeric', 'Age', 'Group']].dropna()
#Spearman correlation because Abeta 42 is a continous variable
spearman_corr_abeta, spearman_p_abeta = stats.spearmanr(abeta_temp_df['Value'], abeta_temp_df['Abeta_42_LCR'])
print(f"Spearman Correlation: ρ = {spearman_corr_abeta:.3f}, p = {spearman_p_abeta:.4f}")
print(f'The sample size is: {len(abeta_temp_df)}')
partial_spearman_abeta = pg.partial_corr(data=abeta_temp_df, x='Value', y='Abeta_42_LCR',
covar='Age', method='spearman')
print("\nSpearman partial correlation:")
print(partial_spearman_abeta)
#Adjust for sex and age
print('='*50)
print('CORRELATION WITH ABETA 42 (ADJUSTED FOR SEX AND AGE)')
print('='*50)
partial_corr_abeta= pg.partial_corr(data= abeta_temp_df, x= 'Abeta_42_LCR', y= 'Value', covar=['Age', 'Sex_numeric'], method= 'spearman')
print(f"ρ: {partial_corr_abeta['p-val'].iloc[0]:.4f}, p = {partial_corr_abeta['p-val'].iloc[0]:.4f}")
print('*'*50)
# Compare unadjusted vs adjusted results
print(f"\nComparison:")
print(f"Unadjusted p-value : {spearman_p_abeta:.4f}")
print(f"Age-adjusted p-value (ANCOVA): {partial_spearman_abeta['p-val'].iloc[0]:.4f}")
if abs(spearman_p_abeta - partial_spearman_abeta['p-val'].iloc[0]) > 0.05:
print("Note: Adjustment for age and sex meaningfully changed the results")
else:
print("Note: Results were robust to age and sex adjustment")
# %%
#Create a scatterplot
plt.figure(figsize=(8,6), dpi= 600)
sns.scatterplot(data= abeta_temp_df, x= "Abeta_42_LCR", y= "Value",
hue= 'Group', palette= 'RdBu', hue_order= group_order)
plt.title(f"hK6 vs. Aβ1-42 (Spearman ρ = {spearman_corr_abeta:.3f}, p = {spearman_p_abeta:.4f})")
plt.xlabel("Aβ1-42 values (pg/mL)")
plt.ylabel("hK6 concentration (ng/mL)")
plt.savefig('comb_hK6vsAbeta42.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correlation with pTau
# %%
#Correlation with p-Tau
print('='*50)
print('CORRELATION WITH P-TAU')
print('='*50)
#Fix dropna issues (the x and y variables are not the same length)
ptau_temp_df = df[["log2_Value","Value", 'P_tau_LCR','Age', 'Sex_numeric', 'Group']].dropna()
#Spearman correlation because Abeta 42 is a continous variable
spearman_corr_ptau, spearman_p_ptau = stats.spearmanr(ptau_temp_df['Value'], ptau_temp_df['P_tau_LCR'])
print(f"Spearman Correlation: ρ = {spearman_corr_ptau:.3f}, p = {spearman_p_ptau:.4f}")
print(f'The sample size is: {len(ptau_temp_df)}')
#Adjust for sex and age
print('='*50)
print('CORRELATION WITH ABETA 42 (ADJUSTED FOR SEX AND AGE)')
print('='*50)
partial_corr_ptau= pg.partial_corr(data= ptau_temp_df, x= 'P_tau_LCR', y= 'Value', covar=['Age', 'Sex_numeric'], method= 'spearman')
print(f"ρ: {partial_corr_ptau['p-val'].iloc[0]:.4f}, p = {partial_corr_ptau['p-val'].iloc[0]:.4f}")
print('*'*50)
# Compare unadjusted vs adjusted results
print(f"\nComparison:")
print(f"Unadjusted p-value : {spearman_p_ptau:.4f}")
print(f"Age-adjusted p-value (ANCOVA): {partial_corr_ptau['p-val'].iloc[0]:.4f}")
if abs(spearman_p_ptau - partial_corr_ptau['p-val'].iloc[0]) > 0.05:
print("Note: Adjustment for age and sex meaningfully changed the results")
else:
print("Note: Results were robust to age and sex adjustment")
print('='*50)
#Create a scatterplot
plt.figure(figsize=(8,6), dpi=600)
sns.scatterplot(data= ptau_temp_df, x= "P_tau_LCR", y= "Value", hue= 'Group',
palette= 'RdBu', hue_order= group_order)
plt.title(f"hK6 vs. p-Tau (Spearman ρ = {spearman_corr_ptau:.3f}, p = {spearman_p_ptau:.4f})")
plt.xlabel("P-tau values (pg/mL)")
plt.ylabel("hK6 concentration (ng/mL)")
plt.savefig('comb_hK6vspTau.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correlation with tTau
# %%
#Correlation with t-Tau
print('='*50)
print('CORRELATION WITH T-TAU')
print('='*50)
#Fix dropna issues (the x and y variables are not the same length)
ttau_temp_df = df[["log2_Value","Value", 'T_tau_LCR','Age', 'Sex_numeric', 'Group']].dropna()
#Spearman correlation because Abeta 42 is a continous variable
spearman_corr_ttau, spearman_p_ttau = stats.spearmanr(ttau_temp_df['Value'], ttau_temp_df['T_tau_LCR'])
print(f"Spearman Correlation: ρ = {spearman_corr_ttau:.3f}, p = {spearman_p_ttau:.4f}")
print(f'The sample size was: {len(ttau_temp_df)}')
partial_spearman_ttau = pg.partial_corr(data=df, x='Value', y='T_tau_LCR',
covar='Age', method='spearman')
print("\nSpearman partial correlation:")
print(partial_spearman_ttau)
#Adjust for sex and age
print('='*50)
print('CORRELATION WITH ABETA 42 (ADJUSTED FOR SEX AND AGE)')
print('='*50)
partial_corr_ttau= pg.partial_corr(data= ttau_temp_df, x= 'T_tau_LCR', y= 'Value', covar=['Age', 'Sex_numeric'], method= 'spearman')
print(f"ρ: {partial_corr_ttau['p-val'].iloc[0]:.4f}, p = {partial_corr_ttau['p-val'].iloc[0]:.4f}")
print('*'*50)
# Compare unadjusted vs adjusted results
print(f"\nComparison:")
print(f"Unadjusted p-value : {spearman_p_ttau:.4f}")
print(f"Age-adjusted p-value (ANCOVA): {partial_corr_ttau['p-val'].iloc[0]:.4f}")
if abs(spearman_p_ttau - partial_corr_ttau['p-val'].iloc[0]) > 0.05:
print("Note: Adjustment for age and sex meaningfully changed the results")
else:
print("Note: Results were robust to age and sex adjustment")
print('='*50)
#Create a scatterplot
plt.figure(figsize=(8,6), dpi= 600)
sns.scatterplot(data= ttau_temp_df, x= "T_tau_LCR", y= "Value",
hue= 'Group', palette= 'RdBu', hue_order= group_order)
plt.title(f"hK6 vs. t-Tau (Spearman ρ = {spearman_corr_ttau:.3f}, p = {spearman_p_ttau:.4f})")
plt.xlabel("T-tau values (pg/mL)")
plt.ylabel("hK6 concentration (ng/mL)")
plt.savefig('comb_hK6vstTau.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correction for multiple testing
# %%
# IMPORTANT
# Correct for multiple testing
print('='*50)
print('MULTIPLE TESTING CORRECTION')
print('='*50)
# Collect all raw p-values
raw_p_values = [
spearman_p_age, partial_spearman_age['p-val'].iloc[0],
p_value_sex, partial_corr_sex['p-val'].iloc[0],
spearman_p_mmse, partial_spearman_mmse['p-val'].iloc[0],
p_value_apoe, ancova_result['p-unc'].iloc[0], # APOE main effect from ANCOVA
spearman_p_abeta, partial_spearman_abeta['p-val'].iloc[0],
spearman_p_ptau, partial_corr_ptau['p-val'].iloc[0],
spearman_p_ttau, partial_corr_ttau['p-val'].iloc[0]
]
# Create a summary DataFrame (FIXED syntax errors)
variables = ['Age', 'Age (adj)',
'Sex', 'Sex (adj)',
'MMSE', 'MMSE (adj)',
'APOE_Status', 'APOE_Status (adj)',
"Aβ1-42", "Aβ1-42 (adj)",
'pTau', 'pTau (adj)',
'tTau', 'tTau (adj)']
tests = ['Spearman', 'Partial Spearman',
'Mann-Whitney U', 'Partial Spearman',
'Spearman', 'Partial Spearman',
'ANOVA/Kruskal-Wallis', 'ANCOVA',
'Spearman', 'Partial Spearman',
'Spearman', 'Partial Spearman',
'Spearman', 'Partial Spearman']
# Apply Holm correction
rejected_Holm, holm_p, _, _ = multipletests(raw_p_values, alpha=0.05, method='holm')
results_summary = pd.DataFrame({
'Variable': variables,
'Test': tests,
'Raw_p_value': raw_p_values,
'Holm_Adjusted_p': holm_p,
'Significant_Holm': rejected_Holm
})
print("Summary of Results with Multiple Testing Correction:")
print(results_summary.round(4))
# Highlight significant results after correction
print("\nSignificant results after Holm correction (p < 0.05):")
significant_results = results_summary[results_summary['Holm_Adjusted_p'] < 0.05]
if len(significant_results) > 0:
print(significant_results[['Variable', 'Test', 'Raw_p_value', 'Holm_Adjusted_p']].round(4))
else:
print("No significant results after multiple testing correction.")
# %% [markdown]
# ## Correlation of hK6 concentration with Age, Sex, MMSE score, APOE status, Abeta, P tau and T tau (non-combined data)
#
# - We used the raw data with Shapiro-Wlilk test, because we are breaking the groups into smaller subgroups, unlike in the earlier analysis. This applies to all analyses herein.
# - MCI_Progressor = MCI_conver_toDem
# %% [markdown]
# ### Correlation with Age
# %%
#Correlation with AGE for MCI non progressors
plot_spearman_correlation('Age', 'Value', data=group_nonprog, xlabel= 'Age (years)',
ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs Age Spearman Correlation for MCI-nonProg")
# %%
#Correlation with AGE for MCI progressors
plot_spearman_correlation('Age', 'Value', data=group_prog, xlabel= 'Age (years)',
ylabel= 'hK6 concentration (ng/mL)', title= "hK6 vs Age Spearman Correlation for MCI-Prog")
# %%
#Correlation with AGE for SCD
plot_spearman_correlation('Age', 'Value', data=group_SCD, xlabel= 'Age (years)',
ylabel= 'hK6 concentration (ng/mL)', title= "hK6 vs Age Spearman Correlation for SCD")
# %%
#Correlation with AGE for AD
plot_spearman_correlation('Age', 'Value', data=group_AD, xlabel= 'Age (years)',
ylabel= 'hK6 concentration (ng/mL)', title= "hK6 vs Age Spearman Correlation for AD")
# %% [markdown]
# ### Correlation with Sex
# %%
#Correlation with Sex MCI-nonProg
print('='*50)
print('CORRELATION WITH SEX FOR MCI-NON PROGRESSORS')
print('='*50)
# Check normality within each group using Shapiro-Wilk
for sex_group in group_nonprog['Sex'].unique():
group_data = group_nonprog[group_nonprog['Sex'] == sex_group]['Value'].dropna()
stat_nonprog, p_nonprog = stats.shapiro(group_data)
print(f"Shapiro-Wilk for {sex_group}: p = {p_nonprog:.4f}")
print(" ")
# %%
print("="*50)
print("The Shapiro-Wilk test rejected the H0, the data deviates from normality")
#Use Mann-Whitney U test, non-normal data
male_data_nonprog = group_nonprog[group_nonprog['Sex']== 'Male']['Value'].dropna()
female_data_nonprog = group_nonprog[group_nonprog['Sex']== 'Female']['Value'].dropna()
stat_sex_nonprog, p_value_sex_nonprog = stats.mannwhitneyu(male_data_nonprog, female_data_nonprog, alternative= 'two-sided')
print(f"\nMann-Whitney U Test: U = {stat_sex_nonprog}, p = {p_value_sex_nonprog:.4f}")
# Create a box plot
plt.figure(figsize=(8, 6), dpi=600)
sns.boxplot(data=group_nonprog, x='Sex', y='Value', hue= "Sex", palette= "RdBu")
sns.stripplot(data=group_nonprog, x='Sex', y='Value', color='black', alpha=0.5, size=4)
plt.xlabel("Sex")
plt.ylabel("hK6 concentration (ng/mL)")
plt.title(f'hK6 Levels by Sex in MCI-nonProg (Mann-Whitney p = {p_value_sex_nonprog:.3f})')
plt.savefig('hK6vsSex_MCI_nonProg.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Correlation with Sex MCI-Prog
print('='*50)
print('CORRELATION WITH SEX FOR MCI-PROGRESSORS')
print('='*50)
# Check normality within each group using Shapiro-Wilk
for sex_group in group_prog['Sex'].unique():
group_data = group_prog[group_prog['Sex'] == sex_group]['Value'].dropna()
stat_prog, p_prog = stats.shapiro(group_data)
print(f"Shapiro-Wilk for {sex_group}: p = {p_prog:.4f}")
# %%
print("The Shapiro-Wilk test rejected the H0, the data deviates from normality")
#Use Mann-Whitney U test, non-normal data
male_data_prog = group_prog[group_prog['Sex']== 'Male']['Value'].dropna()
female_data_prog = group_prog[group_prog['Sex']== 'Female']['Value'].dropna()
stat_sex_prog, p_value_sex_prog = stats.mannwhitneyu(male_data_prog, female_data_prog, alternative= 'two-sided')
print(f"\nMann-Whitney U Test: U = {stat_sex_prog}, p = {p_value_sex_prog:.4f}")
# Create a box plot
plt.figure(figsize=(8, 6), dpi=600)
sns.boxplot(data=group_prog, x='Sex', y='Value', hue= "Sex", palette= "RdBu")
sns.stripplot(data=group_prog, x='Sex', y='Value', color='black', alpha=0.5, size=4)
plt.xlabel("Sex")
plt.ylabel("hK6 concentration (ng/mL)")
plt.title(f'hK6 Levels by Sex in MCI-Prog (Mann-Whitney p = {p_value_sex_prog:.3f})')
plt.savefig('hK6vsSex_MCI_Prog.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Correlation with Sex SCD
print('='*50)
print('CORRELATION WITH SEX FOR SCD')
print('='*50)
# Check normality within each group using Shapiro-Wilk
for sex_group in group_SCD['Sex'].unique():
group_data = group_SCD[group_SCD['Sex'] == sex_group]['Value'].dropna()
stat_SCD, p_SCD = stats.shapiro(group_data)
print(f"Shapiro-Wilk for {sex_group}: p = {p_SCD:.4f}")
# %%
#Use Mann-Whitney U test, non-normal data
male_data_SCD = group_SCD[group_SCD['Sex']== 'Male']['Value'].dropna()
female_data_SCD = group_SCD[group_SCD['Sex']== 'Female']['Value'].dropna()
stat_sex_SCD, p_value_sex_SCD = stats.mannwhitneyu(male_data_SCD, female_data_SCD, alternative= 'two-sided')
print(f"\nMann-Whitney U Test: U = {stat_sex_SCD}, p = {p_value_sex_SCD:.4f}")
# Create a box plot
plt.figure(figsize=(8, 6), dpi=600)
sns.boxplot(data=group_SCD, x='Sex', y='Value', hue= "Sex", palette= "RdBu")
sns.stripplot(data=group_SCD, x='Sex', y='Value', color='black', alpha=0.5, size=4)
plt.xlabel("Sex")
plt.ylabel("hK6 concentration (ng/mL)")
plt.title(f'hK6 Levels by Sex in SCD (Mann-Whitney p = {p_value_sex_SCD:.3f})')
plt.savefig('hK6vsSex_SCD.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Correlation with Sex AD
print('='*50)
print('CORRELATION WITH SEX FOR AD')
print('='*50)
# Check normality within each group using Shapiro-Wilk
for sex_group in group_AD['Sex'].unique():
group_data = group_AD[group_AD['Sex'] == sex_group]['Value'].dropna()
stat_AD, p_AD = stats.shapiro(group_data)
print(f"Shapiro-Wilk for {sex_group}: p = {p_AD:.4f}")
# %%
#Use Mann-Whitney U test, non-normal data
male_data_AD = group_AD[group_AD['Sex']== 'Male']['Value'].dropna()
female_data_AD = group_AD[group_AD['Sex']== 'Female']['Value'].dropna()
stat_sex_AD, p_value_sex_AD = stats.mannwhitneyu(male_data_AD, female_data_AD, alternative= 'two-sided')
print(f"\nMann-Whitney U Test: U = {stat_sex_AD}, p = {p_value_sex_AD:.4f}")
# Create a box plot
plt.figure(figsize=(8, 6), dpi=600)
sns.boxplot(data=group_AD, x='Sex', y='Value', hue= "Sex", palette= "RdBu")
sns.stripplot(data=group_AD, x='Sex', y='Value', color='black', alpha=0.5, size=4)
plt.xlabel("Sex")
plt.ylabel("hK6 concentration (ng/mL)")
plt.title(f'hK6 Levels by Sex in AD (Mann-Whitney p = {p_value_sex_AD:.3f})')
plt.savefig('hK6vsSex_AD.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correlation with MMSE score
# %%
#MCI-nonProgressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
mmse_temp_df_nonprog = group_nonprog[["Value", 'MMSE_PL', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(mmse_temp_df_nonprog)}")
plot_spearman_correlation('MMSE_PL', 'Value', data=mmse_temp_df_nonprog,
xlabel= 'MMSE Score (1-30)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs MMSE Score Spearman Correlation for MCI-nonProg")
# %%
#MCI-progressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
mmse_temp_df_prog = group_prog[["Value", 'MMSE_PL']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(mmse_temp_df_prog)}")
plot_spearman_correlation('MMSE_PL', 'Value', data=mmse_temp_df_prog, xlabel= 'MMSE Score (1-30)',
ylabel= 'hK6 concentration (ng/mL)', title= "hK6 vs MMSE Score Spearman Correlation for MCI-Prog")
# %%
#SCD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
mmse_temp_df_SCD = group_SCD[["Value", 'MMSE_PL']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(mmse_temp_df_SCD)}")
plot_spearman_correlation('MMSE_PL', 'Value', data=mmse_temp_df_SCD, xlabel= 'MMSE Score (1-30)',
ylabel= 'hK6 concentration (ng/mL)', title= "hK6 vs MMSE Score Spearman Correlation for SCD")
# %%
#AD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
mmse_temp_df_AD = group_AD[["Value", 'MMSE_PL']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(mmse_temp_df_AD)}")
plot_spearman_correlation('MMSE_PL', 'Value', data=mmse_temp_df_AD, xlabel= 'MMSE Score (1-30)',
ylabel= 'hK6 concentration (ng/mL)', title= "hK6 vs MMSE Score Spearman Correlation for AD")
# %% [markdown]
# ### Correlation with APOE status
# %%
#Correlation with APOE status for MCI-nonProg
print('='*50)
print("CORRELATION WITH APOE STATUS FOR MCI NON-PROGRESSORS")
print('='*50)
#remove NAs like in the MMSE analysis
apoe_temp_df_nonprog = group_nonprog[["log2_Value","Value", 'APOE', 'Group']].dropna()
# Prepare data for ANOVA - create a list of arrays for each group
apoe_groups = []
for status in apoe_temp_df_nonprog['APOE'].unique():
group_data = apoe_temp_df_nonprog[apoe_temp_df_nonprog['APOE'] == status]['log2_Value']
apoe_groups.append(group_data)
print(f"APOE {status}: n = {len(group_data)}, mean = {np.mean(group_data):.2f} ± {np.std(group_data):.2f}")
# Perform one-way ANOVA
f_stat_apoe_nonprog, p_value_apoe_nonprog = f_oneway(*apoe_groups)
print(f"\nOne-way ANOVA: F = {f_stat_apoe_nonprog:.3f}, p = {p_value_apoe_nonprog:.4f}")
# If ANOVA is significant (p < 0.05), perform Tukey's HSD post-hoc test
if p_value_apoe_nonprog < 0.05:
print("\nPerforming Tukey's HSD post-hoc test:")
# Prepare data for Tukey test
tukey_data = apoe_temp_df_nonprog[['log2_Value', 'APOE']].dropna()
tukey_results = pairwise_tukeyhsd(tukey_data['log2_Value'], tukey_data['APOE'])
print(tukey_results)
# Plot the results
tukey_results.plot_simultaneous()
plt.title("Tukey HSD Post-hoc Test - Confidence Intervals")
plt.show()
else:
print("ANOVA not significant - no post-hoc tests needed.")
# Create a box plot
plt.figure(figsize=(10, 6), dpi= 600)
sns.boxplot(data=apoe_temp_df_nonprog, x='APOE', y='Value', hue = "APOE",palette= 'viridis', boxprops = dict(alpha= 0.6))
sns.stripplot(data=apoe_temp_df_nonprog, x='APOE', y='Value', color='black', alpha=0.5, size=4)
plt.title(f'hK6 Levels by APOE Status in MCI-nonProg (log2 values ANOVA p = {p_value_apoe_nonprog:.3f})')
plt.xticks(rotation=45)
plt.xlabel("APOE Status")
plt.ylabel("hK6 concentration (ng/mL)")
plt.tight_layout()
plt.savefig('hK6vsAPOE_MCI_noprog.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Correlation with APOE status for MCI-Prog
print('='*50)
print("CORRELATION WITH APOE STATUS FOR MCI PROGRESSORS")
print('='*50)
#remove NAs like in the MMSE analysis
apoe_temp_df_prog = group_prog[["log2_Value","Value", 'APOE', 'Group']].dropna()
# Prepare data for ANOVA - create a list of arrays for each group
apoe_groups = []
for status in apoe_temp_df_prog['APOE'].unique():
group_data = apoe_temp_df_prog[apoe_temp_df_prog['APOE'] == status]['log2_Value']
apoe_groups.append(group_data)
print(f"APOE {status}: n = {len(group_data)}, mean = {np.mean(group_data):.2f} ± {np.std(group_data):.2f}")
# Perform one-way ANOVA
f_stat_apoe_prog, p_value_apoe_prog = f_oneway(*apoe_groups)
print(f"\nOne-way ANOVA: F = {f_stat_apoe_prog:.3f}, p = {p_value_apoe_prog:.4f}")
# If ANOVA is significant (p < 0.05), perform Tukey's HSD post-hoc test
if p_value_apoe_prog < 0.05:
print("\nPerforming Tukey's HSD post-hoc test:")
# Prepare data for Tukey test
tukey_data = apoe_temp_df_prog[['log2_Value', 'APOE']].dropna()
tukey_results = pairwise_tukeyhsd(tukey_data['log2_Value'], tukey_data['APOE'])
print(tukey_results)
# Plot the results
tukey_results.plot_simultaneous()
plt.title("Tukey HSD Post-hoc Test - Confidence Intervals")
plt.show()
else:
print("ANOVA not significant - no post-hoc tests needed.")
# Create a box plot
plt.figure(figsize=(10, 6), dpi=600)
sns.boxplot(data=apoe_temp_df_prog, x='APOE', y='Value', hue= "APOE", palette = "viridis", boxprops= dict(alpha=0.6))
sns.stripplot(data=apoe_temp_df_prog, x='APOE', y='Value', color='black', alpha=0.5, size=4)
plt.title(f'hK6 Levels by APOE Status in MCI-Prog (log2 values ANOVA p = {p_value_apoe_prog:.3f})')
plt.xticks(rotation=45)
plt.xlabel("APOE Status")
plt.ylabel("hK6 concentration (ng/mL)")
plt.tight_layout()
plt.savefig('hK6vsAPOE_MCI_prog.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Correlation with APOE status for SCD
print('='*50)
print("CORRELATION WITH APOE STATUS FOR SCD")
print('='*50)
#remove NAs like in the MMSE analysis
apoe_temp_df_SCD = group_SCD[["log2_Value","Value", 'APOE', 'Group']].dropna()
# Prepare data for ANOVA - create a list of arrays for each group
apoe_groups = []
for status in apoe_temp_df_SCD['APOE'].unique():
group_data = apoe_temp_df_SCD[apoe_temp_df_SCD['APOE'] == status]['log2_Value']
apoe_groups.append(group_data)
print(f"APOE {status}: n = {len(group_data)}, mean = {np.mean(group_data):.2f} ± {np.std(group_data):.2f}")
# Perform one-way ANOVA
f_stat_apoe_SCD, p_value_apoe_SCD = f_oneway(*apoe_groups)
print(f"\nOne-way ANOVA: F = {f_stat_apoe_SCD:.3f}, p = {p_value_apoe_SCD:.4f}")
# If ANOVA is significant (p < 0.05), perform Tukey's HSD post-hoc test
if p_value_apoe_SCD < 0.05:
print("\nPerforming Tukey's HSD post-hoc test:")
# Prepare data for Tukey test
tukey_data = apoe_temp_df_SCD[['log2_Value', 'APOE']].dropna()
tukey_results = pairwise_tukeyhsd(tukey_data['log2_Value'], tukey_data['APOE'])
print(tukey_results)
# Plot the results
tukey_results.plot_simultaneous()
plt.title("Tukey HSD Post-hoc Test - Confidence Intervals")
plt.show()
else:
print("ANOVA not significant - no post-hoc tests needed.")
# Create a box plot
plt.figure(figsize=(10, 6), dpi=600)
sns.boxplot(data=apoe_temp_df_SCD, x='APOE', y='Value', hue = "APOE",palette= 'viridis', boxprops = dict(alpha= 0.6))
sns.stripplot(data=apoe_temp_df_SCD, x='APOE', y='Value', color='black', alpha=0.5, size=4)
plt.title(f'hK6 Levels by APOE Status in SCD (log2 values ANOVA p = {p_value_apoe_SCD:.3f})')
plt.xticks(rotation=45)
plt.xlabel("APOE Status")
plt.ylabel("hK6 concentration (ng/mL)")
plt.tight_layout()
plt.savefig('hK6vsAPOE_SCD.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %%
#Correlation with APOE status for AD
print('='*50)
print("CORRELATION WITH APOE STATUS FOR AD")
print('='*50)
#remove NAs like in the MMSE analysis
apoe_temp_df_AD = group_AD[["log2_Value","Value", 'APOE', 'Group']].dropna()
# Prepare data for ANOVA - create a list of arrays for each group
apoe_groups = []
for status in apoe_temp_df_AD['APOE'].unique():
group_data = apoe_temp_df_AD[apoe_temp_df_AD['APOE'] == status]['log2_Value']
apoe_groups.append(group_data)
print(f"APOE {status}: n = {len(group_data)}, mean = {np.mean(group_data):.2f} ± {np.std(group_data):.2f}")
# Perform one-way ANOVA
f_stat_apoe_AD, p_value_apoe_AD = f_oneway(*apoe_groups)
print(f"\nOne-way ANOVA: F = {f_stat_apoe_AD:.3f}, p = {p_value_apoe_AD:.4f}")
# If ANOVA is significant (p < 0.05), perform Tukey's HSD post-hoc test
if p_value_apoe_AD < 0.05:
print("\nPerforming Tukey's HSD post-hoc test:")
# Prepare data for Tukey test
tukey_data = apoe_temp_df_AD[['log2_Value', 'APOE']].dropna()
tukey_results = pairwise_tukeyhsd(tukey_data['log2_Value'], tukey_data['APOE'])
print(tukey_results)
# Plot the results
tukey_results.plot_simultaneous()
plt.title("Tukey HSD Post-hoc Test - Confidence Intervals")
plt.show()
else:
print("ANOVA not significant - no post-hoc tests needed.")
# Create a box plot
plt.figure(figsize=(10, 6), dpi= 600)
sns.boxplot(data=apoe_temp_df_AD, x='APOE', y='Value', hue = "APOE",palette= 'viridis', boxprops = dict(alpha= 0.6))
sns.stripplot(data=apoe_temp_df_AD, x='APOE', y='Value', color='black', alpha=0.5, size=4)
plt.title(f'hK6 Levels by APOE Status in AD (log2 values ANOVA p = {p_value_apoe_AD:.3f})')
plt.xticks(rotation=45)
plt.xlabel("APOE Status")
plt.ylabel("hK6 concentration (ng/mL)")
plt.tight_layout()
plt.savefig('hK6vsAPOE_AD.jpg',
dpi=600,
bbox_inches='tight', # removes extra white space
pad_inches=0.1, # small padding around figure
facecolor='white', # background color
edgecolor='none', # no border
transparent=False) # not transparent
plt.show()
# %% [markdown]
# ### Correlation with Abeta 42
# %%
#MCI-nonProgressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
abeta_temp_df_nonprog = group_nonprog[["Value", 'Abeta_42_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(abeta_temp_df_nonprog)}")
plot_spearman_correlation('Abeta_42_LCR', 'Value', data=abeta_temp_df_nonprog, xlabel= 'Aβ1-42 (pg/mL)',
ylabel= 'hK6 concentration (ng/mL)', title=
"hK6 vs Aβ1-42 Spearman Correlation for MCI-nonProg")
# %%
#MCI-progressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
abeta_temp_df_prog = group_prog[["Value", 'Abeta_42_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(abeta_temp_df_prog)}")
plot_spearman_correlation('Abeta_42_LCR', 'Value', data=abeta_temp_df_prog, xlabel= 'Aβ1-42 (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs Aβ1-42 Spearman Correlation for MCI-Prog")
# %%
#SCD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
abeta_temp_df_SCD = group_SCD[["Value", 'Abeta_42_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(abeta_temp_df_SCD)}")
plot_spearman_correlation('Abeta_42_LCR', 'Value', data=abeta_temp_df_SCD, xlabel= 'Aβ1-42 (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs Aβ1-42 Spearman Correlation for SCD")
# %%
#AD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
abeta_temp_df_AD = group_AD[["Value", 'Abeta_42_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(abeta_temp_df_AD)}")
plot_spearman_correlation('Abeta_42_LCR', 'Value', data=abeta_temp_df_AD, xlabel= 'Aβ1-42 (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs Aβ1-42 Spearman Correlation for AD")
# %% [markdown]
# ### Correlation with p-tau
# %%
#MCI-nonprogressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ptau_temp_df_nonprog = group_nonprog[["Value", 'P_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ptau_temp_df_nonprog)}")
plot_spearman_correlation('P_tau_LCR', 'Value', data=ptau_temp_df_nonprog, xlabel= 'p-Tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)', title=
"hK6 vs p-Tau Spearman Correlation for MCI-nonProg")
# %%
#MCI-progressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ptau_temp_df_prog = group_prog[["Value", 'P_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ptau_temp_df_prog)}")
plot_spearman_correlation('P_tau_LCR', 'Value', data=ptau_temp_df_prog, xlabel= 'p-Tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs p-Tau Spearman Correlation for MCI-Prog")
# %%
#SCD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ptau_temp_df_SCD = group_SCD[["Value", 'P_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ptau_temp_df_SCD)}")
plot_spearman_correlation('P_tau_LCR', 'Value', data=ptau_temp_df_SCD, xlabel= 'p-Tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs p-Tau Spearman Correlation for SCD")
# %%
#AD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ptau_temp_df_AD = group_AD[["Value", 'P_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ptau_temp_df_AD)}")
plot_spearman_correlation('P_tau_LCR', 'Value', data=ptau_temp_df_AD, xlabel= 'p-Tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs p-Tau Spearman Correlation for AD")
# %% [markdown]
# ### Correlation with t-tau
# %%
#MCI-nonprogressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ttau_temp_df_nonprog = group_nonprog[["Value", 'T_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ttau_temp_df_nonprog)}")
plot_spearman_correlation('T_tau_LCR', 'Value', data=ttau_temp_df_nonprog, xlabel= 'T-tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)', title=
"hK6 vs T-tau Spearman Correlation for MCI-nonProg")
# %%
#MCI-progressors
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ttau_temp_df_prog = group_prog[["Value", 'T_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ttau_temp_df_prog)}")
plot_spearman_correlation('T_tau_LCR', 'Value', data=ttau_temp_df_prog, xlabel= 'T-tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs T-tau Spearman Correlation for MCI-Prog")
# %%
#SCD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ttau_temp_df_SCD = group_SCD[["Value", 'T_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ttau_temp_df_SCD)}")
plot_spearman_correlation('T_tau_LCR', 'Value', data=ttau_temp_df_SCD, xlabel= 'T-tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs T-tau Spearman Correlation for SCD")
# %%
#AD
#Because we know there are NAs here, we create a temporary dataframe to avoid errors
ttau_temp_df_AD = group_AD[["Value", 'T_tau_LCR', 'Group']].dropna()
print(f"Original sample size: {len(df)}")
print(f"Sample size after removing missing values: {len(ttau_temp_df_AD)}")
plot_spearman_correlation('T_tau_LCR', 'Value', data=ttau_temp_df_AD, xlabel= 'T-tau (pg/mL)', ylabel= 'hK6 concentration (ng/mL)',
title= "hK6 vs T-tau Spearman Correlation for AD")
# %% [markdown]
# ### Multiple testing correction
# %%
#Make sex binary in the four dataframes
group_prog['Sex_numeric'] = group_prog['Sex'].map({'Male': 0, 'Female': 1}) # or vice versa
group_nonprog['Sex_numeric'] = group_nonprog['Sex'].map({'Male': 0, 'Female': 1}) # or vice versa
group_SCD['Sex_numeric'] = group_SCD['Sex'].map({'Male': 0, 'Female': 1}) # or vice versa
group_AD['Sex_numeric'] = group_AD['Sex'].map({'Male': 0, 'Female': 1}) # or vice versa
# %%
#1. For Age- Spearman correlations
# Perform Spearman correlation (non-parametric, more robust)
spearman_corr_age_nonprog, spearman_p_age_nonprog = stats.spearmanr(group_nonprog['Value'].dropna(), group_nonprog['Age'].dropna())
partial_corr_age_nonprog = pg.partial_corr(data=group_nonprog, x='Age', y='Value', covar=['Sex_numeric'], method='spearman')
spearman_corr_age_prog, spearman_p_age_prog = stats.spearmanr(group_prog['Value'].dropna(), group_prog['Age'].dropna())
partial_corr_age_prog = pg.partial_corr(data=group_prog, x='Age', y='Value', covar=['Sex_numeric'], method='spearman')
spearman_corr_age_SCD, spearman_p_age_SCD = stats.spearmanr(group_SCD['Value'].dropna(), group_SCD['Age'].dropna())
partial_corr_age_SCD = pg.partial_corr(data=group_SCD, x='Age', y='Value', covar=['Sex_numeric'], method='spearman')
spearman_corr_age_AD, spearman_p_age_AD = stats.spearmanr(group_AD["Value"].dropna(), group_AD['Age'].dropna())
partial_corr_age_AD = pg.partial_corr(data=group_AD, x='Age', y='Value', covar=['Sex_numeric'], method='spearman')
# %%
#2. For Sex- Mann-Whitney U tests, ALREADY HAVE IT
# %%
#3. For MMSE score- spearman correlations
group_nonprog_mmse = group_nonprog.dropna(subset= ['MMSE_PL'])
group_prog_mmse = group_prog.dropna(subset= ['MMSE_PL'])
group_SCD_mmse = group_SCD.dropna(subset= ['MMSE_PL'])
group_AD_mmse = group_AD.dropna(subset= ['MMSE_PL'])
spearman_corr_MMSE_nonprog, spearman_p_MMSE_nonprog = stats.spearmanr(group_nonprog_mmse['Value'].dropna(), group_nonprog_mmse['MMSE_PL'].dropna())
partial_corr_MMSE_nonprog = pg.partial_corr(data=group_nonprog_mmse, x='MMSE_PL', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_MMSE_prog, spearman_p_MMSE_prog = stats.spearmanr(group_nonprog_mmse['Value'].dropna(), group_nonprog_mmse['MMSE_PL'].dropna())
partial_corr_MMSE_prog = pg.partial_corr(data=group_prog_mmse, x='MMSE_PL', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_MMSE_SCD, spearman_p_MMSE_SCD = stats.spearmanr(group_nonprog_mmse['Value'].dropna(), group_nonprog_mmse['MMSE_PL'].dropna())
partial_corr_MMSE_SCD = pg.partial_corr(data=group_SCD_mmse, x='MMSE_PL', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_MMSE_AD, spearman_p_MMSE_AD = stats.spearmanr(group_nonprog_mmse["Value"].dropna(), group_nonprog_mmse['MMSE_PL'].dropna())
partial_corr_MMSE_AD = pg.partial_corr(data=group_AD_mmse, x='MMSE_PL', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
# %%
#4. For APOE status- ANOVA, ALREADY HAVE IT
# %%
#5. For Abeta 42- Spearman correlations
group_nonprog_abeta = group_nonprog.dropna(subset= ['Abeta_42_LCR'])
group_prog_abeta = group_prog.dropna(subset= ['Abeta_42_LCR'])
group_SCD_abeta = group_SCD.dropna(subset= ['Abeta_42_LCR'])
group_AD_abeta = group_AD.dropna(subset= ['Abeta_42_LCR'])
spearman_corr_Abeta_nonprog, spearman_p_Abeta_nonprog = stats.spearmanr(group_nonprog_abeta['Value'].dropna(), group_nonprog_abeta['Abeta_42_LCR'].dropna())
partial_corr_abeta_nonprog = pg.partial_corr(data=group_nonprog_abeta, x='Abeta_42_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_Abeta_prog, spearman_p_Abeta_prog = stats.spearmanr(group_prog_abeta['Value'].dropna(), group_prog_abeta['Abeta_42_LCR'].dropna())
partial_corr_abeta_prog = pg.partial_corr(data=group_prog_abeta, x='Abeta_42_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_Abeta_SCD, spearman_p_Abeta_SCD = stats.spearmanr(group_SCD_abeta['Value'].dropna(), group_SCD_abeta['Abeta_42_LCR'].dropna())
partial_corr_abeta_SCD = pg.partial_corr(data=group_SCD_abeta, x='Abeta_42_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_Abeta_AD, spearman_p_Abeta_AD = stats.spearmanr(group_AD_abeta["Value"].dropna(), group_AD_abeta['Abeta_42_LCR'].dropna())
partial_corr_abeta_AD = pg.partial_corr(data=group_AD_abeta, x='Abeta_42_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
# %%
#6. For p-tau- Spearman correlations
group_nonprog_ptau = group_nonprog.dropna(subset= ['P_tau_LCR'])
group_prog_ptau = group_prog.dropna(subset= ['P_tau_LCR'])
group_SCD_ptau = group_SCD.dropna(subset= ['P_tau_LCR'])
group_AD_ptau = group_AD.dropna(subset= ['P_tau_LCR'])
spearman_corr_ptau_nonprog, spearman_p_ptau_nonprog = stats.spearmanr(group_nonprog_ptau['Value'].dropna(), group_nonprog_ptau['P_tau_LCR'].dropna())
partial_corr_ptau_nonprog = pg.partial_corr(data=group_nonprog_ptau, x='P_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_ptau_prog, spearman_p_ptau_prog = stats.spearmanr(group_prog_ptau['Value'].dropna(), group_prog_ptau['P_tau_LCR'].dropna())
partial_corr_ptau_prog = pg.partial_corr(data=group_prog_ptau, x='P_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_ptau_SCD, spearman_p_ptau_SCD = stats.spearmanr(group_SCD_ptau['Value'].dropna(), group_SCD_ptau['P_tau_LCR'].dropna())
partial_corr_ptau_SCD = pg.partial_corr(data=group_SCD_ptau, x='P_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_ptau_AD, spearman_p_ptau_AD = stats.spearmanr(group_AD_ptau["Value"].dropna(), group_AD_ptau['P_tau_LCR'].dropna())
partial_corr_ptau_AD = pg.partial_corr(data=group_AD_ptau, x='P_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
# %%
#7. For t-tau- Spearman correlations
group_nonprog_ttau = group_nonprog.dropna(subset= ['T_tau_LCR'])
group_prog_ttau = group_prog.dropna(subset= ['T_tau_LCR'])
group_SCD_ttau = group_SCD.dropna(subset= ['T_tau_LCR'])
group_AD_ttau = group_AD.dropna(subset= ['T_tau_LCR'])
spearman_corr_ttau_nonprog, spearman_p_ttau_nonprog = stats.spearmanr(group_nonprog_ttau['Value'].dropna(), group_nonprog_ttau['T_tau_LCR'].dropna())
partial_corr_ttau_nonprog = pg.partial_corr(data=group_nonprog_ttau, x='T_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_ttau_prog, spearman_p_ttau_prog = stats.spearmanr(group_prog_ttau['Value'].dropna(), group_prog_ttau['T_tau_LCR'].dropna())
partial_corr_ttau_prog = pg.partial_corr(data=group_prog_ttau, x='T_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_ttau_SCD, spearman_p_ttau_SCD = stats.spearmanr(group_SCD_ttau['Value'].dropna(), group_SCD_ttau['T_tau_LCR'].dropna())
partial_corr_ttau_SCD = pg.partial_corr(data=group_SCD_ttau, x='T_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
spearman_corr_ttau_AD, spearman_p_ttau_AD = stats.spearmanr(group_AD_ttau["Value"].dropna(), group_AD_ttau['T_tau_LCR'].dropna())
partial_corr_ttau_AD = pg.partial_corr(data=group_AD_ttau, x='T_tau_LCR', y='Value', covar=['Sex_numeric', 'Age'], method='spearman')
# %%
print('='*50)
print('MULTIPLE TESTING CORRECTION FOR NON-COMBINED RESULTS')
print('='*50)
raw_p_values_noncomb = [spearman_p_age_nonprog, spearman_p_age_prog, spearman_p_age_SCD, spearman_p_age_AD,
partial_corr_age_nonprog['p-val'].iloc[0], partial_corr_age_prog['p-val'].iloc[0], partial_corr_age_SCD['p-val'].iloc[0], partial_corr_age_AD['p-val'].iloc[0],
p_value_sex_nonprog, p_value_sex_prog, p_value_sex_SCD, p_value_sex_AD,
spearman_p_MMSE_nonprog, spearman_p_MMSE_prog, spearman_p_MMSE_SCD, spearman_p_MMSE_AD,
partial_corr_MMSE_nonprog['p-val'].iloc[0], partial_corr_MMSE_prog['p-val'].iloc[0], partial_corr_MMSE_SCD['p-val'].iloc[0], partial_corr_MMSE_AD['p-val'].iloc[0],
p_value_apoe_nonprog, p_value_apoe_prog, p_value_apoe_SCD, p_value_apoe_AD,
spearman_p_Abeta_nonprog, spearman_p_Abeta_prog, spearman_p_Abeta_SCD, spearman_p_Abeta_AD,
partial_corr_abeta_nonprog['p-val'].iloc[0], partial_corr_abeta_prog['p-val'].iloc[0], partial_corr_abeta_SCD['p-val'].iloc[0], partial_corr_abeta_AD['p-val'].iloc[0],
spearman_p_ptau_nonprog, spearman_p_ptau_prog, spearman_p_ptau_SCD, spearman_p_ptau_AD,
spearman_p_ttau_nonprog, spearman_p_ttau_prog, spearman_p_ttau_SCD, spearman_p_ttau_AD ] # [Age, Sex, MMSE, APOE, Abeta, pTau, tTau]
#Apply Holm correction
rejected, holm_p_noncomb,_,_ = multipletests(raw_p_values_noncomb, alpha= 0.05, method= 'holm')
# Create a summary DataFrame
variables_noncomb = ['Age','Age', 'Age', 'Age'
,'Age (adj)','Age (adj)', 'Age (adj)', 'Age (adj)'
,'Sex','Sex', 'Sex', 'Sex'
,'MMSE','MMSE', 'MMSE', 'MMSE'
,'MMSE (adj)','MMSE (adj)', 'MMSE (adj)', 'MMSE (adj)'
,'APOE_Status', 'APOE_Status', 'APOE_Status', 'APOE_Status'
,"Aβ1-42", "Aβ1-42", "Aβ1-42", "Aβ1-42"
,"Aβ1-42 (adj)", "Aβ1-42 (adj)", "Aβ1-42 (adj)", "Aβ1-42 (adj)"
,'pTau', 'pTau', 'pTau', 'pTau'
,'tTau', 'tTau', 'tTau', 'tTau']
tests_noncomb = ['Spearman','Spearman','Spearman','Spearman',
'Adjusted Spearman','Adjusted Spearman','Adjusted Spearman','Adjusted Spearman',
'Mann-Whitney U','Mann-Whitney U','Mann-Whitney U','Mann-Whitney U',
'Spearman','Spearman','Spearman','Spearman',
'Adjusted Spearman','Adjusted Spearman','Adjusted Spearman','Adjusted Spearman',
'ANOVA','ANOVA','ANOVA','ANOVA',
'Spearman','Spearman','Spearman','Spearman',
'Adjusted Spearman','Adjusted Spearman','Adjusted Spearman','Adjusted Spearman',
'Spearman','Spearman','Spearman','Spearman',
'Spearman','Spearman','Spearman','Spearman']
groups_noncomb = ['MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD',
'MCI non Progressor', 'MCI Progressor', 'SCD', 'AD']
results_summary_noncomb = pd.DataFrame({
'Variable': variables_noncomb,
'Test': tests_noncomb,
'Group': groups_noncomb,
'Raw_p_value': raw_p_values_noncomb,
'Holm_p': holm_p_noncomb
})
print("Summary of Results with Multiple Testing Correction:")
print(results_summary_noncomb.round(4))
# Highlight significant results after correction
print("\nSignificant results after Holm correction (p < 0.05):")
significant_results_noncomb = results_summary_noncomb[results_summary_noncomb['Holm_p'] < 0.05]
if len(significant_results_noncomb) > 0:
print(significant_results_noncomb[['Group','Variable', 'Holm_p']])
else:
print("No significant results after multiple testing correction.")
# %% [markdown]
# ## Kruskal Wallis for all variables
# - Since we do not have to test for equal variances/normality, it is easier to perform Kruskal-Wallis on these variables.
# - MMSE Score variables has to be analyzed with Kruskal-Wallis test, due to it being an ordinal variable.
# %%
# Rename column 'Value' to 'hK6' in-place
df.rename(columns={'Value': 'hK6'}, inplace=True)
# %%
variables = ['Age', 'MMSE_PL', 'Abeta_42_LCR', 'P_tau_LCR', 'T_tau_LCR', 'hK6'] #Variables to perform Kruskal Wallis to
# Create a results dataframe to store all statistics
results = []
print("Kruskal-Wallis Test Results for Multiple Variables")
print("=" * 60)
for var in variables:
print(f"\nVariable: {var}")
print("-" * 30)
# Extract data for each group
group_data = []
for group in group_order:
group_values = df[df['Group'] == group][var].dropna()
group_data.append(group_values)
print(f"Group {group}: n={len(group_values)}, median={group_values.median():.2f}")
# Perform Kruskal-Wallis test
h_stat, p_value = stats.kruskal(*group_data)
print(f"Kruskal-Wallis H-statistic: {h_stat:.3f}, p-value: {p_value:.4f}")
# Store results
result_dict = {
'Variable': var,
'H_statistic': h_stat,
'p_value': p_value,
'significant': p_value < 0.05
}
# Add group medians and sample sizes
for i, group in enumerate(group_order):
result_dict[f'Group_{group}_median'] = group_data[i].median()
result_dict[f'Group_{group}_n'] = len(group_data[i])
results.append(result_dict)
# Perform Mann-Whitney U post-hoc tests with Holm correction if significant
if p_value < 0.05:
print("Significant difference found! Performing Mann-Whitney U post-hoc tests:")
# Calculate number of comparisons for correction
n_groups = len(group_order)
n_comparisons = n_groups * (n_groups - 1) // 2
print(f"Performing {n_comparisons} pairwise comparisons (Holm correction)")
print("\nPairwise comparisons (Mann-Whitney U):")
print("-" * 50)
posthoc_results = []
raw_p_values_for_Holm = [] # store raw p-values for Holm correction
# First pass: collect all raw p-values
for i in range(n_groups):
for j in range(i + 1, n_groups):
group1, group2 = group_order[i], group_order[j]
data1 = group_data[i]
data2 = group_data[j]
# Perform Mann-Whitney U test
stat, p_raw = stats.mannwhitneyu(data1, data2, alternative='two-sided')
# Calculate effect size (rank-biserial correlation)
n1, n2 = len(data1), len(data2)
effect_size = 1 - (2 * stat) / (n1 * n2) # simple effect size measure
result = {
'Comparison': f'{group1} vs {group2}',
'U_statistic': stat,
'Raw_p_value': p_raw,
'Effect_size': effect_size
}
posthoc_results.append(result)
raw_p_values_for_Holm.append(p_raw)
# Apply Holm correction to all p-values (AFTER collecting them all)
rejected, holm_p_values, _, _ = multipletests(raw_p_values_for_Holm, alpha=0.05, method='holm')
# Update results with Holm-corrected p-values
for k, result in enumerate(posthoc_results):
result['Adjusted_p_value'] = holm_p_values[k]
result['Significant'] = holm_p_values[k] < 0.05
result['Comparison_index'] = k
# Calculate significance marker first to avoid f-string issues
sig_marker = '*' if result['Significant'] else ''
# Use consistent quotes and correct variable names
print(f"{result['Comparison']}: U = {result['U_statistic']:.1f}, "
f"raw p = {result['Raw_p_value']:.4f}, "
f"adj p = {result['Adjusted_p_value']:.4f} {sig_marker}")
# Add posthoc results to the main results
result_dict['posthoc_results'] = pd.DataFrame(posthoc_results)
# Convert results to dataframe
results_df = pd.DataFrame(results)
# Display summary table
print("\n" + "=" * 80)
print("SUMMARY TABLE")
print("=" * 80)
summary_cols = ['Variable', 'H_statistic', 'p_value', 'significant']
for group in group_order:
summary_cols.extend([f'Group_{group}_median', f'Group_{group}_n'])
print(results_df[summary_cols].round(4))
# Create visualizations
# use correct group order
n_vars = len(variables)
n_cols = min(3, n_vars)
n_rows = (n_vars + n_cols - 1) // n_cols # ceiling division
fig, axes = plt.subplots(n_rows, n_cols, figsize=(5*n_cols, 5*n_rows), dpi=600)
if n_vars > 1:
axes = axes.flatten()
else:
axes = [axes]
for i, var in enumerate(variables):
if i < len(axes):
# Create boxplot (for the subplots)
sns.boxplot(data=df, x='Group', y=var, ax=axes[i], hue='Group',
palette='RdBu', legend=False, hue_order=group_order, order=group_order)
sns.stripplot(data=df, x='Group', y=var, ax=axes[i], legend=False,
alpha=0.6, color='black', size=3, order=group_order)
p_val = results_df[results_df['Variable'] == var]['p_value'].values[0]
sig_stars = '***' if p_val < 0.001 else '**' if p_val < 0.01 else '*' if p_val < 0.05 else 'ns'
axes[i].set_title(f'{var}\np = {p_val:.4f} {sig_stars}', size=13)
axes[i].set_xlabel('Groups', size=12)
axes[i].set_ylabel(var)
axes[i].set_xticks(range(len(df['Group'].unique())))
axes[i].set_xticklabels(group_order, ha="right", rotation_mode="anchor")
axes[i].tick_params(axis="x", rotation=45)
# Hide any empty subplots
for j in range(len(variables), len(axes)):
axes[j].set_visible(False)
plt.tight_layout()
plt.suptitle('Distribution of Variables Across the Four Diagnostic Groups', fontsize=16, y=1.02)
plt.savefig('kruskal_wallis_results_7_figures.jpg', dpi=600, bbox_inches='tight')
plt.show()
# Save results to CSV
results_df.to_csv('kruskal_wallis_results.csv', index=False)
print("\nResults saved to 'kruskal_wallis_results.csv'")
# Print significant findings
print("\n" + "=" * 80)
print("SIGNIFICANT FINDINGS")
print("=" * 80)
for _, row in results_df.iterrows():
if row['significant']:
print(f"\n{row['Variable']} shows significant differences across groups (p = {row['p_value']:.4f})")
# FIXED: Check if posthoc_results exists and is not empty
if 'posthoc_results' in row and hasattr(row['posthoc_results'], 'empty'):
posthoc_df = row['posthoc_results']
if not posthoc_df.empty:
# Make sure the column name matches what you stored
sig_column = 'Significant' if 'Significant' in posthoc_df.columns else 'significant'
sig_comparisons = posthoc_df[posthoc_df[sig_column] == True]
if not sig_comparisons.empty:
print("Significant pairwise differences (Holm-corrected):")
for _, comp in sig_comparisons.iterrows():
print(f" {comp['Comparison']}: adj p = {comp['Adjusted_p_value']:.4f}")
else:
print(" No significant pairwise differences found after Holm correction")
else:
print(" No posthoc results available (empty DataFrame)")
else:
print(" No posthoc results available")
# 1. ANCOVA for continuous variables (age-adjusted group comparisons)
age_sensitive_vars = ['MMSE_PL', 'Abeta_42_LCR', 'P_tau_LCR', 'T_tau_LCR', 'hK6']
print("AGE-ADJUSTED COMPARISONS (ANCOVA)")
print("=" * 50)
for var in age_sensitive_vars:
print(f"\nVariable: {var} (adjusted for age)")
print("-" * 30)
# Perform ANCOVA with age as covariate
ancova_result = pg.ancova(data=df, dv=var, between='Group', covar='Age')
print(ancova_result)
# Compare with Kruskal-Wallis result
kw_p = results_df[results_df['Variable'] == var]['p_value'].values[0]
ancova_p = ancova_result['p-unc'].iloc[0] # Group effect p-value
print(f"Kruskal-Wallis p: {kw_p:.4f}")
print(f"ANCOVA (age-adjusted) p: {ancova_p:.4f}")
if abs(kw_p - ancova_p) > 0.05:
print("→ Age appears to be an important confounder")
else:
print("→ Results are robust to age adjustment")
# 2. Alternatively, use partial correlations for continuous relationships
print("\n" + "=" * 50)
print("AGE-ADJUSTED CORRELATIONS WITH GROUP (Partial Correlation)")
print("=" * 50)
# Create a numeric group variable for correlation
df['Group_numeric'] = df['Group'].map({group: i for i, group in enumerate(group_order)})
for var in age_sensitive_vars:
partial_corr = pg.partial_corr(data=df, x='Group_numeric', y=var,
covar=['Age'], method='spearman')
print(f"{var}: ρ = {partial_corr['r'].iloc[0]:.3f}, p = {partial_corr['p-val'].iloc[0]:.4f}")
# %% [markdown]
# ## Extra section: creation of tables for patient statistics, not necessary
# %%
#define for InterQuartile Range
def calculate_iqr(dataframe):
q1 = dataframe.quantile(0.25)
q3 = dataframe.quantile(0.75)
return q3 - q1
def q1(dataframe):
return dataframe.quantile(0.25)
def q3(dataframe):
return dataframe.quantile(0.75)
# %%
#Interquartile range
iqr = df.groupby("Group")["Value"].apply(calculate_iqr)
print(iqr)
print("-"*30)
iqr_log2 = df.groupby("Group")["log2_Value"].apply(calculate_iqr)
print(iqr_log2)
# %%
#Find the Q1 and Q3 by applying the functions
q_25 = df.groupby("Group")["Value"].apply(q1) #25th quantile
print(q_25)
print("-"*30)
q_75 = df.groupby("Group")['Value'].apply(q3) #75th quantile
print(q_75)
print("-"*30)
q_25_log = df.groupby("Group")["log2_Value"].apply(q1)
print(q_25_log)
print("-"*30)
q_75_log = df.groupby("Group")['log2_Value'].apply(q3)
print(q_75_log)
# %%
# Group by 'Group' and apply multiple aggregation functions to 'Value' and 'log2_Value
# Save the tables to an Excel file
table1_values = df.groupby('Group').agg({'Value': ['mean', 'median', 'std', 'min', 'max']})
Table1 = pd.concat([table1_values,q_25, q_75], join= "inner", axis=1) #concat to add quartiles
Table1.columns = ["Mean", "Median", "std", "Min", "Max", "Q1", "Q3"]#add headers to q_25 and q_75
print(Table1)
Table1.to_excel("my_table_values.xlsx")
table1_values_log2 = df.groupby('Group').agg({'log2_Value': ['mean', 'median', 'std', 'min', 'max']})
Table1_log2 = pd.concat([table1_values_log2,q_25_log, q_75_log], join= "inner", axis=1) #concat
Table1_log2.columns = ["Mean", "Median", "std", "Min", "Max", "Q1", "Q3"]#add headers to q_25 and q_75
print(Table1_log2)
Table1_log2.to_excel("my_table_log2_values.xlsx")

KLK6_ELISA_data_analysis.ipynb at commit eb12641, under MIT · at the source

Overview

Authors: Miyo K Chatanaka¹, Antoninus Soosaipillai², Amanda Cano^3,4, Adelina Orellana^3,4, Mercè Boada^3,4, Ioannis Prassas^1,5, Xavier Morató^3,4, Eleftherios P Diamandis²

Institutions: University of Toronto (Canada); Mount Sinai Hospital (Canada); Lunenfeld-Tanenbaum Research Institute (Canada); Instituto de Salud Carlos III (Spain); Universitat Internacional de Catalunya (Spain); Biomedical Research Networking Center on Neurodegenerative Diseases (Spain); University Health Network (Canada)

Journal: Clinical proteomics, volume 23, issue 1, article 23

Dates: received 26 September 2025; accepted 3 December 2025; published online 15 March 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1186/s12014-025-09577-x · PMID 41834054 · PMCID PMC13104293 · OpenAlex W7136088233

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), Alzheimer's / dementia (population), clinical / translational (subfield)

Methods: Statistics, Connectivity

Keywords: Human kallikrein 6, Alzheimer’s disease, Dementia biomarkers, Amyloid beta 1–42, Total tau, Phosphorylated tau

Topic: Coagulation, Bradykinin, Polyphosphates, and Angioedema (Genetics, Medicine), according to OpenAlex

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

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

miyohtnk/PhD_KLK6_data_analysis

License: MIT

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Commit: eb12641861466c13c8629529eb7f29084b7fbbe8, 26 September 2025

Languages: Jupyter (1)

Size: 5 files, 1 script

Software Heritage: not archived

Found in: the text, “Statistical analysis”

Holds: README, license file, 1 notebook

Not found: 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 30 September 2026: the link answers

30 September 2026: the link answers

1 file

KLK6_ELISA_data_analysis.ipynb, Jupyter, 2,452 lines, 2 matches
repository limit reached (2,000 files or 30 MB): the rest is at the source (2 files)

At the source: github.com/miyohtnk/PhD_KLK6_data_analysis

miyohtnk/phd

License: none: the authors keep all their rights

State: the link is dead, verified on 30 September 2026

Evidence: found in the paper

Software Heritage: not archived

Found in: “Data availability”

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

Availability: 1 check, the latest on 30 September 2026: the link is dead

30 September 2026: the link is dead

At the source: github.com/miyohtnk/PhD/_KLK6/_data/_analysis

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;
1 script, each with its path and the digest of its content;
2 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.

Code and data availability statement

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

it points to the authors' code: miyohtnk/phd

Read it in the paper: doi.org/10.1186/s12014-025-09577-x.

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 6 keywords, 50 references.

Cite

This paper

Chatanaka, M. K., Soosaipillai, A., Cano, A., Orellana, A., Boada, M., Prassas, I., Morató, X., & Diamandis, E. P. (2026). Validation of human kallikrein 6 in the cerebrospinal fluid of patients with progressive and non-progressive alzheimer's disease: correlation with other biomarkers. Clinical proteomics, 23(1), 23. https://doi.org/10.1186/s12014-025-09577-x

BibTeX

@article{chatanaka2026validation, author = {Chatanaka, Miyo K and Soosaipillai, Antoninus and Cano, Amanda and Orellana, Adelina and Boada, Mercè and Prassas, Ioannis and Morató, Xavier and Diamandis, Eleftherios P}, title = {{Validation of human kallikrein 6 in the cerebrospinal fluid of patients with progressive and non-progressive alzheimer's disease: correlation with other biomarkers}}, journal = {Clinical proteomics}, year = {2026}, month = mar, volume = {23}, number = {1}, pages = {23}, publisher = {BMC}, issn = {1542-6416}, doi = {10.1186/s12014-025-09577-x}, url = {https://doi.org/10.1186/s12014-025-09577-x}, pmid = {41834054}, pmcid = {PMC13104293} }

RIS

TY - JOUR AU - Chatanaka, Miyo K AU - Soosaipillai, Antoninus AU - Cano, Amanda AU - Orellana, Adelina AU - Boada, Mercè AU - Prassas, Ioannis AU - Morató, Xavier AU - Diamandis, Eleftherios P TI - Validation of human kallikrein 6 in the cerebrospinal fluid of patients with progressive and non-progressive alzheimer's disease: correlation with other biomarkers T2 - Clinical proteomics J2 - Clin Proteomics PY - 2026 DA - 2026/03/15 VL - 23 IS - 1 SP - 23 SN - 1542-6416 PB - BMC DO - 10.1186/s12014-025-09577-x UR - https://doi.org/10.1186/s12014-025-09577-x LA - en ER -

CSL-JSON

{ "id": "10.1186/s12014-025-09577-x", "type": "article-journal", "title": "Validation of human kallikrein 6 in the cerebrospinal fluid of patients with progressive and non-progressive alzheimer's disease: correlation with other biomarkers", "container-title": "Clinical proteomics", "author": [ { "family": "Chatanaka", "given": "Miyo K" }, { "family": "Soosaipillai", "given": "Antoninus" }, { "family": "Cano", "given": "Amanda" }, { "family": "Orellana", "given": "Adelina" }, { "family": "Boada", "given": "Mercè" }, { "family": "Prassas", "given": "Ioannis" }, { "family": "Morató", "given": "Xavier" }, { "family": "Diamandis", "given": "Eleftherios P" } ], "container-title-short": "Clin Proteomics", "volume": "23", "issue": "1", "page": "23", "DOI": "10.1186/s12014-025-09577-x", "PMID": "41834054", "PMCID": "PMC13104293", "ISSN": "1542-6416", "publisher": "BMC", "URL": "https://doi.org/10.1186/s12014-025-09577-x", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 15 ] ] } }

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

Similar papers

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

[1] doi:10.1038/s41398-026-04081-8 [code]: Functional system-specific brain aging across the Alzheimer's disease continuum.
Journal: Translational psychiatry
In common: Pingouin, statsmodels, seaborn, 5 other tools, Alzheimer's / dementia, clinical / translational, 1 reference
[2] doi:10.1371/journal.pone.0343722 [code]: Comprehensive methodology for sample enrichment in EEG biomarker studies for Alzheimer's risk classification.
Journal: PloS one
In common: Pingouin, statsmodels, seaborn, 5 other tools, Alzheimer's / dementia, clinical / translational
[3] doi:10.64898/2026.05.06.26352540 [code]: Generating synthetic tau-PET scans in Alzheimer’s disease from MRI, blood biomarkers and demographics with deep learning
Journal: medRxiv (preprint)
In common: seaborn, scikit-learn, pandas, 3 other tools, Alzheimer's / dementia, clinical / translational, 2 references
[4] doi:10.1002/alz.71649 [code]: Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: Pingouin, statsmodels, seaborn, 5 other tools, Alzheimer's / dementia
[5] doi:10.1002/alz.71365 [code]: Benchmarking speech biomarkers of Alzheimer's against cognitive and neural measures.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: Pingouin, statsmodels, seaborn, 5 other tools, Alzheimer's / dementia
[6] doi:10.34133/csbj.0042 [code]: Using Steady-State Visual Evoked Potentials to Characterize Wide-Ranging Retinopathy Linked to <i>CRB1</i>: Implications for Clinical Trials.
Journal: Computational and structural biotechnology journal
In common: Pingouin, statsmodels, seaborn, 5 other tools, clinical / translational
[7] doi:10.1038/s43587-026-01096-0 [code]: Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis.
Journal: Nature aging
In common: Pingouin, statsmodels, seaborn, 5 other tools, Alzheimer's / dementia
[8] doi:10.1162/imag.a.1259 [code]: Neuronal avalanches as a predictive biomarker for guiding tailored BCI training programs.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Pingouin, statsmodels, seaborn, 5 other tools, clinical / translational
[9] doi:10.1038/s41531-026-01380-1 [code]: Identifying maximal beta power from directional subthalamic local field potentials in Parkinson's disease.
Journal: NPJ Parkinson's disease
In common: Pingouin, statsmodels, seaborn, 5 other tools, clinical / translational
[10] doi:10.7554/elife.108673 [code]: Adaptive behavior is guided by integrated representations of controlled and non-controlled information.
Journal: eLife
In common: Pingouin, statsmodels, seaborn, 5 other tools

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.