Subgingival microbiota composition is associated with brain health in the general population-the PAROMIND study.

Code ↔ Paper

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

[1] § Methods › Topological data analysis framework › SAFE: enrichment analysis ↔ analysis.ipynb, lines 29–92 · score 0.99 · Geriatric Depression Scale, diastolic blood pressure, Mini Mental State, blood triglycerides, MEDAS score, MIND score
[2] § Methods › Phenotyping › Cognitive and mental health assessments ↔ analysis.ipynb, lines 29–92 · score 0.82 · Geriatric Depression Scale, Mini Mental State, Recall, Vocabulary, Word, inverted
[3] § Methods › Topological data analysis framework › SAFE: enrichment analysis ↔ analysis.ipynb, lines 596–662 · score 0.81 · robust Aitchison distance, forward model selection, ordiR2step, matrix, vegan, variance
[4] § Methods › Group analysis ↔ analysis.ipynb, lines 1533–1548 · score 0.75 · diastolic blood pressure, BMI, HDL, HbA1c, cholesterol, triglycerides

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,057 lines · 68 KB · no license · 4 matches

# %% [markdown]
# # Import libraries and data
# %%
import numpy as np
np.int = np.int32
from sklearn.preprocessing import MinMaxScaler
from sklearn.cluster import DBSCAN
from tmap.tda import mapper, Filter
from tmap.tda.cover import Cover
from tmap.tda.metric import Metric
from tmap.tda.utils import optimize_dbscan_eps
import pandas as pd
import networkx as nx
# %%
from pathlib import Path
code_dir=Path.cwd()
project_dir=code_dir.parent
input_dir=project_dir/"input"
output_dir=project_dir/"output/tda/"
tmp_dir=project_dir/"tmp"
output_dir.mkdir(exist_ok=True, parents=True)
# %% [markdown]
# # Prepare data
#
# %%
# Grouping variables in categories
demographics = [
"demographics_age",
"demographics_sex",
"demographics_education_isced",
]
cognitive_scores = [
"cognition_g_factor_inverted",
"cognition_tmt_a_inverted",
"cognition_tmt_b_inverted",
'cognition_vocabulary_b_test',
'cognition_word_list_recall',
'cognition_animal_naming_test',
'cognition_mini_mental_state_exam',
]
neuropsychiatric_scores = [
'psych_phq9_sum',
'psych_phq15_sum',
'psych_geriatrics_depression_scale',
'psych_gad7_sum',
]
cardiovascular_risk_factors = [
'cvrisk_systolic_blood_pressure_mmhg',
'cvrisk_diastolic_blood_pressure_mmhg',
'cvrisk_BMI',
'cvrisk_smoking_currently',
"blood_cholesterol_mg_dl",
"blood_hdl_mg_dl",
"blood_ldl_mg_dl",
"blood_triglycerides_mg_dl",
"blood_hba1c",
]
inflammation = [
"inflammation_hsCRP",
"inflammation_leukocytes",
]
paro_variables = [
"paro_cal_mean",
"paro_dmft",
"paro_plaqueindex",
"paro_bop"
]
diet_scores = [
"diet_mind_score",
"diet_medas_score",
"diet_dash1_score"
]
microbiome_variables = oral_microbiome_genus.columns.tolist()
imaging_means = ['imaging_thickness_volume_mean','imaging_fw_gm_mean', 'imaging_fw_wm_mean', 'imaging_fat_gm_mean', 'imaging_fat_wm_mean']
metadata_variables = demographics + cognitive_scores + neuropsychiatric_scores + cardiovascular_risk_factors + inflammation + paro_variables + imaging_means + diet_scores
# %%
# Setup naming dictionary (shortened for brevity)
variable_styling_dict = {
"base_age":"Age",
"...":"...",
}
# %%
# Load taxa abundance data, sample metadata
metadata = pd.read_csv(input_dir/"data/metadata_df.csv", index_col=0)
oral_microbiome_genus = pd.read_csv(input_dir/"data/oral_microbiome_genus.csv", index_col=0)
X = oral_microbiome_genus
metadata = metadata.loc[metadata.index.isin(X.index)][metadata_variables]
X = X.loc[X.index.isin(metadata.index)]
# %%
# z-score non-microbiome data
from sklearn.preprocessing import StandardScaler
def is_binary(series):
return set(series.dropna().unique()) <= {0, 1}
binary_columns = [col for col in metadata.columns if is_binary(metadata[col])]
metadata[binary_columns] = metadata[binary_columns].astype('object')
continuous_columns = [col for col in metadata.columns if not is_binary(metadata[col])]
metadata[continuous_columns] = StandardScaler().fit_transform(metadata[continuous_columns])
# %%
metadata_categories = [col.split("_")[0] for col in metadata.columns.tolist()]
microbiome_categories = ["genus"] * len(X.columns.to_list())
# %% [markdown]
# # Aitchison distance
# %%
from biom.table import Table
data = X.T.to_numpy()
samples = X.index.to_list()
observation = X.columns.to_list()
table = Table(data, observation, samples)
from gemelli.rpca import rpca
_, distance = rpca(table)
dm = distance.to_data_frame()
# %% [markdown]
# # Mapper
# %%
# Initiate a Mapper object
tm = mapper.Mapper(verbose=1)
# Projection
metric = Metric(metric="precomputed")
lens = [Filter.PCOA(components=[0, 1], metric=metric, random_state=100)]
projected_X = tm.filter(dm, lens=lens)
# Covering, clustering & mapping
eps = optimize_dbscan_eps(X, threshold=95)
clusterer = DBSCAN(eps=eps, min_samples=3)
cover = Cover(projected_data=MinMaxScaler().fit_transform(projected_X), resolution=30, overlap=1.5)
graph = tm.map(data=X, cover=cover, clusterer=clusterer)
print(graph.info())
# %%
# Apply spring embedding
initial_nodepos = {idx:graph.nodePos[idx] for idx in range(graph.nodePos.shape[0])}
pos = nx.spring_layout(graph, k = 0.2, pos = initial_nodepos, seed=42)
graph.nodePos = np.array([pos[key] for key in pos.keys()])
for idx, node in enumerate(graph.nodes):
graph.nodes[idx]["pos"] = pos[idx].tolist()
import pickle
pickle.dump(graph, open(output_dir/"graph.pickle", "wb"))
edgelist_3col = nx.to_pandas_edgelist(graph)
edgelist_3col["dist"] = 1
edgelist_3col.to_csv(output_dir/"mapper_graph_3col.txt", sep="\t", index=False, header=None)
# %% [markdown]
# # SAFE
# %%
# Setup function to transform data from subject to node level
def transform2node_data(graph, data, mode='mean'):
map_fun = {'sum': np.sum,
"mean": np.nanmean}
if mode not in ["sum", "mean"]:
raise SyntaxError('Wrong provided parameters.')
else:
aggregated_fun = map_fun[mode]
nodes = graph.nodes
dv = data.values
if data is not None:
node_data = {nid: aggregated_fun(dv[attr['sample'], :], 0)
for nid, attr in nodes.items()}
node_data = pd.DataFrame.from_dict(node_data,
orient='index',
columns=data.columns)
return node_data
# %%
node_subject_mapping = {node:list(graph.nodes[idx]["sample_names"]) for idx,node in enumerate(graph.nodes)}
def transform2node_data_bin(df, node_subject_mapping):
# Computes node level percentage of binary variables
t = [(k, x) for k, v in node_subject_mapping.items() for x in v]
hierarchical_index = pd.MultiIndex.from_tuples(t)
df_long = pd.DataFrame(index=hierarchical_index, columns=df.columns)
df_long.index.rename(["node","subject"], inplace=True)
for node,sub in df_long.index:
df_long.loc[(node,sub),:] = df.loc[sub,:]
for col in df_long.columns :
df_long[col] = df_long[col].astype(object)
df_long.reset_index(inplace=True)
df_transformed = df_long[["node"] + df.columns.tolist()].groupby("node").agg(np.sum) / df_long[["node"] + df.columns.tolist()].groupby("node").agg(len)
return df_transformed
# %%
# Prepare data to be used by safepy
metadata_transformed = transform2node_data(graph, metadata.drop(['demographics_sex', 'cvrisk_smoking_currently'], axis=1), mode="mean")
binary_transformed = transform2node_data_bin(metadata_plus_imaging[['demographics_sex', 'cvrisk_smoking_currently']], node_subject_mapping)
metadata_transformed = metadata_transformed.join(binary_transformed)
oral_microbiome_genus_transformed = transform2node_data(graph, oral_microbiome_genus, mode="mean")
data_transformed = metadata_transformed.join(oral_microbiome_genus_transformed)
data_transformed.to_csv(output_dir/"mapper_graph_metadata.txt", sep="\t", index=True)
# %%
# Perform SAFE
from safepy import safe
sf = safe.SAFE(path_to_safe_data=f"{output_dir}/safe/")
sf.random_seed = 0
sf.load_network(network_file=f"{output_dir}/mapper_graph_3col.txt")
sf.load_attributes(attribute_file=f"{output_dir}/mapper_graph_metadata.txt")
sf.define_neighborhoods()
num_permutations = 5000
sf.compute_pvalues(num_permutations=num_permutations)
# %%
network_enrichment_scores = pd.DataFrame(sf.nes, columns=data_transformed.columns)
network_enrichment_scores_signif = pd.DataFrame(sf.nes_binary, columns=data_transformed.columns)
network_enrichment_scores_signif_pos = (network_enrichment_scores > 0) & (network_enrichment_scores_signif)
network_enrichment_scores_signif_neg = (network_enrichment_scores < 0) & (network_enrichment_scores_signif)
# %%
nes_min = network_enrichment_scores.min().min()
nes_max = network_enrichment_scores.max().max()
print(nes_min, nes_max)
# %%
safe_summary = sf.attributes.copy()
safe_summary.drop("id", axis=1)
safe_summary.set_index("name", inplace=True)
# %%
# Compute enrichment ratio + extract positive and negative enrichment
# enrichment ratio = n of enriched neighborhoods / n of nodes
safe_summary["enrichment_ratio"] = safe_summary["num_neighborhoods_enriched"] / len(graph.nodes)
safe_summary["num_neighborhoods_enriched_pos"] = pd.DataFrame(network_enrichment_scores_signif_pos.sum(axis=0), index = safe_summary.index)
safe_summary["num_neighborhoods_enriched_neg"] = pd.DataFrame(network_enrichment_scores_signif_neg.sum(axis=0), index = safe_summary.index)
safe_summary["category"] = [idx[0] for idx in safe_summary.index.str.split("_")]
safe_summary.to_csv(output_dir/'metadata_safe_summary.csv')
# %% [markdown]
# # Draw PCA based on SAFE score
# %%
color_codes = {
'genus': '#68D391', # A fresh green
'cognition': '#D53F8C', # A vibrant pink
'psych': '#FFD700', # Gold
'cvrisk': '#1E90FF', # DodgerBlue
'blood': '#1E90FF', # dark red
'inflammation': '#9E2A2B', # dark red
"diet": "#036c5f", #Teal
'paro': '#DD6B20', # A distinctive orange
'demographics': '#808080', # DarkSlateGray
'imaging': '#8A2BE2', # BlueViolet
}
color_codes_light = {
'genus': '#a3ebc7', # pastel green
'cognition': '#f086b1', # pastel pink
'psych': '#ffeb99', # pastel gold/yellow
'cvrisk': '#80c4ff', # pastel blue
'blood': '#80c4ff', # pastel blue
'inflammation': '#d16b6b', # pastel red
"diet": "#4d9e8f", #pastel teal
'paro': '#ffae73', # pastel orange
'demographics': '#b3b3b3', # pastel gray
'imaging': '#c393f7' # pastel violet
}
# %%
import plotly.offline as py_offline
import plotly.graph_objs as go
from sklearn.decomposition import PCA
from sklearn.preprocessing import MinMaxScaler
import pandas as pd
# Perform PCA
pca = PCA()
pca_result = pca.fit_transform(network_enrichment_scores.T)
# Select top enrichments for metadata and oral microbiome genus
top_10_col_sorter = lambda df, col: df.sort_values(col, ascending=False).head(10).index
top_metadata = top_10_col_sorter(safe_summary.loc[data_transformed.columns], 'enrichment_ratio')
top_genus = top_10_col_sorter(safe_summary.loc[oral_microbiome_genus.columns], 'enrichment_ratio')
# Normalize enrichment ratio
mx_scale = MinMaxScaler(feature_range=(10, 30))
scaled_enrichment = mx_scale.fit_transform(safe_summary[['enrichment_ratio']])
# Prepare data for plotting
data = []
categories = safe_summary['category'].unique()
for cat in categories:
mask = safe_summary['category'] == cat
category_data = safe_summary[mask]
scaled_size = scaled_enrichment[mask].flatten()
trace = go.Scatter(
x=pca_result[mask, 0],
y=pca_result[mask, 1],
mode='markers',
name=variable_styling_dict[cat],
marker=dict(
color=color_codes[cat],
size=10,#scaled_size, #
opacity=0.5
),
text=[variable_styling_dict[cat] if not "genus" in cat else " ".join(cat.split("_")[1:]) for cat in category_data.index],
textfont=dict(
size=12 # specify the size of the text
),
)
data.append(trace)
# Layout configuration
layout = go.Layout(
xaxis=dict(title=f"PC1 ({pca.explained_variance_ratio_[0] * 100:.2f}%)"),
yaxis=dict(title=f"PC2 ({pca.explained_variance_ratio_[1] * 100:.2f}%)"),
width=1200, height=500,
title="",
font=dict(size=15),
hovermode='closest'
)
# Plot using offline mode in Plotly
fig = go.Figure(data=data, layout=layout)
py_offline.plot(fig, auto_open=False)
fig.write_html(output_dir/"enrichment_analysis/PCA.html")
fig.write_image(output_dir/"enrichment_analysis/PCA.png", format="png", scale=10)
fig.write_image(output_dir/"enrichment_analysis/PCA.svg", format="svg")
pyo.iplot(fig, config={'responsive': True}) # Use iframe as the renderer
# %%
import plotly.offline as py_offline
import plotly.graph_objs as go
from sklearn.decomposition import PCA
from sklearn.preprocessing import MinMaxScaler
import pandas as pd
import numpy as np
def avoid_text_overlap(x, y, text_list, offset=0.02):
new_positions = []
seen_positions = set()
for i in range(len(x)):
new_x, new_y = x[i], y[i]
while (new_x, new_y) in seen_positions:
new_y += offset
seen_positions.add((new_x, new_y))
new_positions.append((new_x, new_y))
return new_positions
# Perform PCA
pca = PCA()
pca_result = pca.fit_transform(network_enrichment_scores.T)
# Select top enrichments for metadata and oral microbiome genus
top_10_col_sorter = lambda df, col: df.sort_values(col, ascending=False).head(10).index
top_metadata = top_10_col_sorter(safe_summary.loc[data_transformed.columns], 'enrichment_ratio')
top_genus = top_10_col_sorter(safe_summary.loc[oral_microbiome_genus.columns], 'enrichment_ratio')
# Normalize enrichment ratio
mx_scale = MinMaxScaler(feature_range=(10, 30))
scaled_enrichment = mx_scale.fit_transform(safe_summary[['enrichment_ratio']])
# Prepare data for plotting
data = []
categories = safe_summary['category'].unique()
for cat in categories:
mask = safe_summary['category'] == cat
category_data = safe_summary[mask]
scaled_size = scaled_enrichment[mask].flatten()
x_coords = pca_result[mask, 0]
y_coords = pca_result[mask, 1]
text_labels = [
variable_styling_dict[cat] if not "genus" in cat else " ".join(cat.split("_")[1:])
for cat in category_data.index
]
# Adjust text positions to avoid overlap
adjusted_positions = avoid_text_overlap(x_coords, y_coords, text_labels)
trace = go.Scatter(
x=[pos[0] for pos in adjusted_positions],
y=[pos[1] for pos in adjusted_positions],
mode='markers+text',
name=variable_styling_dict[cat],
marker=dict(
color=color_codes[cat],
size=10, # scaled_size,
opacity=0.5
),
text=text_labels,
textfont=dict(
size=10 # specify the size of the text
),
showlegend = False
)
data.append(trace)
# Layout configuration
layout = go.Layout(
xaxis=dict(title=f"PC1 ({pca.explained_variance_ratio_[0] * 100:.2f}%)"),
yaxis=dict(title=f"PC2 ({pca.explained_variance_ratio_[1] * 100:.2f}%)"),
width=1322*1.5, height=794*1.5,
title="",
font=dict(size=15),
hovermode='closest',
template="plotly_white"
)
# Plot using offline mode in Plotly
fig = go.Figure(data=data, layout=layout)
py_offline.plot(fig, auto_open=False)
# Uncomment next line if interactive plot within a notebook is required
py_offline.iplot(fig)
fig.write_html(str(output_dir / "enrichment_analysis/PCA_annotated.html"))
fig.write_image(str(output_dir / "enrichment_analysis/PCA_annotated.png"), format="png", scale=10)
fig.write_image(str(output_dir / "enrichment_analysis/PCA_annotated.svg"), format="svg")
# %% [markdown]
# # envfit
# %%
Path(output_dir/"envfit/").mkdir(exist_ok=True)
X.to_csv(output_dir/'envfit/microbiome.csv',index=True)
metadata.astype(float).to_csv(output_dir/'envfit/metadata.csv',index=True)
# %%
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
import numpy as np
from statsmodels.sandbox.stats.multicomp import multipletests
importr("vegan")
def envfit_metadata(genus_path,metadata_path):
rcode = """
genus_table <- read.csv('{path_data}',row.names = 1)
metadata <- read.csv('{path_metadata}',row.names = 1,check.names=FALSE)
dist <- vegdist(genus_table, method="robust.aitchison")
ord <- capscale(dist ~ -1)
""".format(path_data=genus_path,path_metadata=metadata_path)
robjects.r(rcode)
envfit_result = robjects.r(
"""
fit <- envfit(ord,metadata,permutations = 5000, na.rm=TRUE)
fit$vectors
""")
fit_result = pd.DataFrame(columns=["r2","pvals","Source","End"],index=envfit_result[envfit_result.names.index("arrows")].rownames)
fit_result.loc[:,"r2"] = envfit_result[envfit_result.names.index("r")]
fit_result.loc[:, "pvals"] = envfit_result[envfit_result.names.index("pvals")]
fit_result.loc[:, ["Source","End"]] = np.array(envfit_result[envfit_result.names.index("arrows")])
return fit_result
import time
t1 = time.time()
envfit_df = envfit_metadata(f'{output_dir}/envfit/microbiome.csv',
f'{output_dir}/envfit/metadata.csv',
)
envfit_df["pvals_fdr"] = multipletests(envfit_df["pvals"], method='fdr_bh')[1]
print('envfit takes', time.time() - t1)
# %%
n_envfit = X.shape[0]
p_envfit = 1
envfit_df["adjusted_r2"] = 1 - (1 - envfit_df["r2"]) * ((n_envfit - 1) / (n_envfit - p_envfit -1))
# %% [markdown]
# # adonis
# %%
Path(output_dir/"adonis/").mkdir(exist_ok=True)
X.to_csv(output_dir/'adonis/microbiome.csv',index=True)
# %%
# Perform PERMANOVA with adonis
import rpy2.robjects.pandas2ri as rpypandas
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
importr("vegan")
def run_adonis(genus_path, metadata_path):
"""
Perform cumulative PERMANOVA using adonis to assess the variance explained by all metadata covariates.
Returns a pandas DataFrame with adonis results.
"""
r_code = f"""
genus_table <- read.csv('{genus_path}', row.names = 1)
metadata <- read.csv('{metadata_path}', row.names = 1, check.names = FALSE)
dist <- vegdist(genus_table, method = "robust.aitchison")
adonis_result <- adonis2(dist ~ ., data = metadata, na.action = na.omit, permutations = 5000,
by = "margin"
)
list(adonis_table = adonis_result)
"""
r_output = robjects.r(r_code)
adonis_table_r = r_output.rx2("adonis_table")
adonis_results_df = rpypandas.rpy2py_dataframe(adonis_table_r)
return adonis_results_df
adonis_df = pd.DataFrame(columns=["r2","pvals"])
for col in metadata.columns:
metadata[col].to_csv(output_dir/'adonis/metadata.csv',index=True)
fit_result = run_adonis(f'{output_dir}/adonis/microbiome.csv', f'{output_dir}/adonis/metadata.csv')
adonis_df.loc[col,"r2"] = fit_result.loc[col,"R2"]
adonis_df.loc[col,"pvals"] = fit_result.loc[col,"Pr(>F)"]
adonis_df.loc["pvals_fdr"] = multipletests(adonis_df["pvals"], method='fdr_bh')[1]
adonis_df.to_csv(output_dir/'adonis/adonis_results.csv', index=True)
# %%
n_adonis = X.shape[0]
p_adonis = 1
adonis_df["adjusted_r2"] = 1 - (1 - adonis_df["r2"]) * ((n_adonis - 1) / (n_adonis - p_adonis -1))
# %%
# Merge enrichment ratio, envfit and adonis results into a single table
compared_table = safe_summary.loc[metadata.columns].sort_values('enrichment_ratio',ascending=False)
neg_log_p_thresh = -np.log10(0.05)
compared_table.loc[:,'envfit_adjusted_r2'] = list(envfit_df.loc[compared_table.index,'adjusted_r2'])
compared_table.loc[:,'envfit_p_fdr'] = list(envfit_df.loc[compared_table.index,'pvals_fdr'])
compared_table.loc[:,'envfit_neg_log_p_fdr'] = -np.log10(compared_table.loc[:,'envfit_p_fdr'])
compared_table.loc[:, 'envfit_signif'] = compared_table['envfit_p_fdr'] < 0.05
compared_table.loc[:,'adonis_adjusted_r2'] = list(adonis_df.loc[compared_table.index,'adjusted_r2'])
compared_table.loc[:,'adonis_p_fdr'] = list(adonis_df.loc[compared_table.index,'pvals_fdr'])
compared_table.loc[:,'adonis_neg_log_p_fdr'] = -np.log10(compared_table.loc[:,'adonis_p_fdr'])
compared_table.loc[:, 'adonis_signif'] = compared_table['adonis_p_fdr'] < 0.05
compared_table = compared_table.fillna(0)
compared_table.index = [variable_styling_dict[var] for var in compared_table.index]
compared_table.to_csv(output_dir/'enrichment_analysis/comparison_table.csv',index=True)
# %%
from scipy.stats import spearmanr
print(spearmanr(compared_table["enrichment_ratio"], compared_table["envfit_adjusted_r2"]))
print(spearmanr(compared_table["enrichment_ratio"], compared_table["adonis_adjusted_r2"]))
# %% [markdown]
# # ordiR2step
# %%
import pandas as pd
import rpy2.robjects.pandas2ri as rpypandas
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
from pathlib import Path
# Ensure the R package 'vegan' is imported
importr("vegan")
def run_ordiR2step(genus_path, metadata_path, output_path):
"""
Performs forward model selection using ordiR2step to find the best set of predictors.
This function calculates the robust Aitchison distance and then uses ordiR2step
to build a cumulative model, adding variables one by one based on their
contribution to the explained variance (R-squared).
The final table of the best model is saved to a CSV file.
Args:
genus_path (str): Path to the microbiome genus data CSV file.
metadata_path (str): Path to the full metadata CSV file with all potential predictors.
output_path (str): Path to save the resulting ANOVA table CSV.
"""
print("Performing forward model selection with ordiR2step...")
# R code to be executed
r_code = """
# 1. Load data
genus_table <- read.csv('{genus_path}', row.names = 1)
metadata <- read.csv('{metadata_path}', row.names = 1, check.names=FALSE)
# 2. Calculate the distance matrix
dist <- vegdist(genus_table, method="robust.aitchison")
# 3. Define the null and full models for the selection scope
# Null model (starting point)
mod0 <- capscale(dist ~ 1, data=metadata, na.action=na.omit)
# Full model
mod1 <- capscale(dist ~ ., data=metadata, na.action=na.omit)
# 4. Perform the forward selection using ordiR2step
step_result <- ordiR2step(mod0, scope = formula(mod1), direction = "forward", permutations = 5000)
# 5. Extract the ANOVA table from the final selected model
final_anova <- as.data.frame(step_result$anova)
# Return the final table
final_anova
""".format(genus_path=genus_path, metadata_path=metadata_path)
# Execute the R code
final_model_r = robjects.r(r_code)
# Convert the R dataframe to a pandas dataframe
final_model_df = rpypandas.rpy2py_dataframe(final_model_r)
# Save the results
final_model_df.to_csv(output_path, index=True)
print(f"ordiR2step results saved to {output_path}")
return final_model_df
# %%
from sklearn.impute import KNNImputer # We impute missing data as ordiR2step does not allow NaNs
ordiR2step_dir = Path("./output/ordir2step")
ordiR2step_dir.mkdir(exist_ok=True)
imp = KNNImputer()
imp.fit(metadata)
metadata_imp = imp.transform(metadata)
metadata_imp = pd.DataFrame(metadata_imp, index=metadata.index, columns=metadata.columns)
metadata_imp[analysis_columns].to_csv(ordiR2step_dir/'metadata.csv',index=True)
X[X.index.isin(metadata_imp.index)].to_csv(ordiR2step_dir/'microbiome.csv',index=True)
genus_data_path = str(ordiR2step_dir / 'microbiome.csv')
full_metadata_path = str(ordiR2step_dir / 'metadata.csv')
results_output_path = str(ordiR2step_dir / 'ordir2step_results.csv')
ordir2_model_results = run_ordiR2step(genus_data_path, full_metadata_path, results_output_path)
# %% [markdown]
# # Plot enrichment ratios and adonis results
# %%
from plotly import tools
import plotly.graph_objs as go
fig = tools.make_subplots(rows=5, cols=1, shared_xaxes=True, vertical_spacing=0.07, subplot_titles=['SAFE', 'envfit', '', 'adonis', ''])
plotting_table = compared_table.copy()
plotting_table.index = [idx.split(" ")[0] for idx in plotting_table.index]
plotting_table = plotting_table.sort_values(by=['enrichment_ratio'], ascending=False)
# Plotting enrichment ratio
fig.append_trace(
go.Bar(
y=plotting_table.loc[:, 'enrichment_ratio'] * 100,
x=plotting_table.index,
marker=dict(
color=[color_codes[plotting_table.loc[fea, 'category']] for fea in plotting_table.index],
line=dict(width=1)
),
orientation='v',
showlegend=False
), 1, 1
)
# Plotting envfit_adjusted_r2
fig.append_trace(
go.Bar(
y=plotting_table.loc[:, 'envfit_adjusted_r2'],
x=plotting_table.index,
marker=dict(
color=[
color_codes[plotting_table.loc[fea, 'category']] if plotting_table.loc[fea, 'envfit_p_fdr'] < 0.05 else color_codes_light[plotting_table.loc[fea, 'category']]
for fea in plotting_table.index
],
line=dict(width=1)
),
orientation='v',
showlegend=False
), 2, 1
)
# Plotting envfit_neg_log_p_fdr
fig.append_trace(
go.Bar(
y=plotting_table.loc[:, 'envfit_neg_log_p_fdr'],
x=plotting_table.index,
marker=dict(
color=[
color_codes[plotting_table.loc[fea, 'category']] if plotting_table.loc[fea, 'envfit_p_fdr'] < 0.05 else color_codes_light[plotting_table.loc[fea, 'category']]
for fea in plotting_table.index
],
line=dict(width=1)
),
orientation='v',
showlegend=False
), 3, 1
)
# Add a line for threshold
fig.add_shape(
type="line",
x0=-0.5, y0=-np.log10(0.05), x1=len(plotting_table.index), y1=-np.log10(0.05),
line=dict(color="gray", width=1, dash="dash"),
xref='x1', yref='y3'
)
# Plotting adonis_adjusted_r2
fig.append_trace(
go.Bar(
y=plotting_table.loc[:, 'adonis_adjusted_r2'],
x=plotting_table.index,
marker=dict(
color=[
color_codes[plotting_table.loc[fea, 'category']] if plotting_table.loc[fea, 'adonis_p_fdr'] < 0.05 else color_codes_light[plotting_table.loc[fea, 'category']]
for fea in plotting_table.index
],
line=dict(width=1)
),
orientation='v',
showlegend=False
), 4, 1
)
# Plotting adonis_neg_log_p_fdr
fig.append_trace(
go.Bar(
y=plotting_table.loc[:, 'adonis_neg_log_p_fdr'],
x=plotting_table.index,
marker=dict(
color=[
color_codes[plotting_table.loc[fea, 'category']] if plotting_table.loc[fea, 'adonis_p_fdr'] < 0.05 else color_codes_light[plotting_table.loc[fea, 'category']]
for fea in plotting_table.index
],
line=dict(width=1)
),
orientation='v',
showlegend=False
), 5, 1
)
# Add a line for threshold
fig.add_shape(
type="line",
x0=-0.5, y0=-np.log10(0.05), x1=len(plotting_table.index), y1=-np.log10(0.05),
line=dict(color="gray", width=1, dash="dash"),
xref='x1', yref='y5'
)
# Adjust layout
fig.layout.xaxis1.title = ""
fig.layout.yaxis1.title = "Enrichment (%)"
fig.layout.yaxis2.title = "R2adj"
fig.layout.yaxis3.title = "-log10(pFDR)"
fig.layout.yaxis4.title = "R2adj"
fig.layout.yaxis5.title = "-log10(pFDR)"
fig.layout.margin.t = 40
fig.layout.height = 600
fig.layout.width = 1000
fig.update_xaxes(tickangle=45)
# Save the plots
fig.write_html(output_dir / "enrichment_analysis/barplot_enrichment_metadata_vertical3.html")
fig.write_image(output_dir / "enrichment_analysis/barplot_enrichment_metadata_vertical3.png", format="png", scale=10)
fig.write_image(output_dir / "enrichment_analysis/barplot_enrichment_metadata_vertical3.svg", format="svg")
fig.show()
# %%
# Plot barplot of enrichment ratios only
from plotly import tools
import plotly.graph_objs as go
fig = tools.make_subplots(rows=1, cols=1, shared_xaxes=True, vertical_spacing=0.15, subplot_titles=[''])
plotting_table = compared_table.copy()
plotting_table.index = [idx.split(" ")[0] for idx in plotting_table.index]
plotting_table = plotting_table.sort_values(by=['enrichment_ratio'], ascending=False)
# Plotting enrichment ratio
fig.append_trace(
go.Bar(
y=plotting_table.loc[:, 'enrichment_ratio'] * 100,
x=plotting_table.index,
marker=dict(
color=[color_codes[plotting_table.loc[fea, 'category']] for fea in plotting_table.index],
line=dict(width=1)
),
orientation='v',
showlegend=False
), 1, 1
)
# Adjust layout
fig.layout.xaxis1.title = ""
fig.layout.yaxis1.title = "Enrichment (%)"
fig.layout.margin.t = 40
fig.layout.height = 350
fig.layout.width = 1000
fig.update_xaxes(tickangle=45)
# Save the plots
fig.write_html(output_dir / "enrichment_analysis/barplot_enrichment_metadata_vertical.html")
fig.write_image(output_dir / "enrichment_analysis/barplot_enrichment_metadata_vertical.png", format="png", scale=10)
fig.write_image(output_dir / "enrichment_analysis/barplot_enrichment_metadata_vertical.svg", format="svg")
# Show the plots
fig.show()
# Save the resulting table
compared_table.to_csv(output_dir / "enrichment_analysis/compared_bar_result_data.csv")
# %% [markdown]
# # Plot ordiR2step results
# %%
ordir2_plot_table = ordir2_model_results.copy()
ordir2_plot_table = ordir2_plot_table.rename(index=variable_styling_dict_plus)
individual_r2 = ordir2_plot_table['R2.adj'].diff().fillna(ordir2_plot_table['R2.adj'])
variables = ordir2_plot_table.index
fig = go.Figure()
for i, var in enumerate(variables):
legend_entry = f"{var} (+{individual_r2.iloc[i]*100:.2f}%)"
fig.add_trace(go.Bar(
y=[''],
x=[individual_r2.iloc[i]],
name=legend_entry,
orientation='h'
))
fig.update_layout(
barmode='stack',
title_text='',
xaxis_title="Adjusted R²",
yaxis_title="",
legend_title_text="Significant factors",
legend=dict(traceorder='normal')
)
fig.write_html(output_dir / "enrichment_analysis/stacked_barplot_ordir2_legend_values.html")
fig.write_image(output_dir / "enrichment_analysis/stacked_barplot_ordir2_legend_values.png", format="png", scale=10)
fig.write_image(output_dir / "enrichment_analysis/stacked_barplot_ordir2_legend_values.svg", format="svg")
fig.show()
# %% [markdown]
# # Plotting enrichment landscapes
# %%
import plotly.graph_objs as go
def plot_enrichment_landscape(
graph=None, # networkx graph object | list of node position (e.g., [0.69,0.77]) as node attribute 'pos'.
attribute=None,
network_enrichment_scores_signif=None,
fade_nonsignificant_nodes=True,
variable=None, # str, variable to plot
node_colormap="balance", # colormap of node coloring
title=None, # title to display on the plot
titlefont_size=15, # size of the title font
annotation_text="", # text to display as annotation in the plot
color_range_min=-5, # range of colorbar
color_range_max=5, # range of colorbar
nonsignif_opacity=0.4, # opacity of nonsignificant nodes
node_line_width = 0.7, # width of lines around nodes
node_line_color = "darkgray", # color of lines around nodes, can also be list of length n_nodes to color each node line individually
show_colorbar=True, # indicate whether to display the colorbar
colorbar_annotation_text="", # annotation text of the colorbar
width=500, # figure width
height=500, # figure height
):
G = graph.copy()
assert len(attribute) == len(G.nodes), "len(attribute) does not equal len(graph.nodes())"
signif_idx = network_enrichment_scores_signif[network_enrichment_scores_signif[variable] == 1].index.tolist()
if fade_nonsignificant_nodes == True:
opacity = [1 if idx in signif_idx else nonsignif_opacity for idx in range(len(graph.nodes))]
else:
opacity = 1
node_line_color_list = node_line_color
node_text = []
for idx, _ in enumerate(G.nodes):
node_text.append(f'Node: {idx} Value: {attribute[idx]:.3f} Subjects: {", ".join(G.nodes[_]["sample_names"].tolist())}')
edge_x = []
edge_y = []
for edge in G.edges():
x0, y0 = G.nodes[edge[0]]['pos']
x1, y1 = G.nodes[edge[1]]['pos']
edge_x.append(x0)
edge_x.append(x1)
edge_x.append(None)
edge_y.append(y0)
edge_y.append(y1)
edge_y.append(None)
edge_trace = go.Scatter(
x=edge_x, y=edge_y,
line=dict(width=0.5, color='#888'),
hoverinfo='none',
mode='lines')
node_x = []
node_y = []
for node in G.nodes():
x, y = G.nodes[node]['pos']
node_x.append(x)
node_y.append(y)
node_trace = go.Scatter(
x=node_x, y=node_y,
mode='markers',
hoverinfo='text',
marker=dict(
showscale=show_colorbar,
colorscale=node_colormap,
reversescale=True,
color=[],
cmin=color_range_min,
cmax=color_range_max,
size=10,
colorbar=dict(
thickness=15,
title=f'{colorbar_annotation_text}',
xanchor='left',
titleside='right',
),
line_width=2))
node_trace.marker.color = attribute
node_trace.marker.opacity = opacity
node_trace.text = node_text
node_trace.marker.line["width"] = node_line_width
node_trace.marker.line["color"] = node_line_color_list
fig = go.Figure(data=[edge_trace, node_trace],
layout=go.Layout(
width=width, height=height,
title={
"text":f'{title}',
"x":0.5,
"y":0.95,
},
titlefont_size=titlefont_size,
showlegend=False,
hovermode='closest',
plot_bgcolor='white', # Background color for the plotting area
paper_bgcolor='white', # Background color for the entire figure
margin=dict(b=20,l=5,r=5,t=40),
annotations=[ dict(
text=f"{annotation_text}",
showarrow=False,
xref="paper", yref="paper",
x=0.005, y=-0.002 ) ],
xaxis=dict(showgrid=False, zeroline=False, showticklabels=False),
yaxis=dict(showgrid=False, zeroline=False, showticklabels=False))
)
return fig
# %%
# Plot network without annotation
from plotly.offline import plot
import plotly.io as pio
attribute = [1] * metadata_transformed.shape[0]
fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=False,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable="demographics_age", # pass as dummy variable so that code works
node_colormap="Blues_r",
title="Microbiome network",
titlefont_size=40,
annotation_text="",
color_range_min=None,
color_range_max=None,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=False,
colorbar_annotation_text="",
height=500,
width=450,
)
fig.show()
pio.write_image(fig, output_dir/f"enrichment_analysis/network.png")
plot(fig, filename=Path(output_dir/f"enrichment_analysis/network.html").as_posix(), auto_open=False)
fig.write_image(output_dir/f"enrichment_analysis/network.svg")
# %%
# Plot enrichment landscapes of non-microbiome phenotypes
for variable in list(metadata.columns):
attribute = network_enrichment_scores[variable]
fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=True,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable=variable,
node_colormap="balance",
title=variable_styling_dict[variable].split(" ")[0],
titlefont_size=25,
annotation_text="",
color_range_min=-10,
color_range_max=10,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=False,
colorbar_annotation_text="",
height=500,
width=450,
)
pio.write_image(fig, output_dir/f"enrichment_landscapes/network_{variable}.png", scale=10)
pio.write_image(fig, output_dir/f"enrichment_landscapes/network_{variable}.pdf", scale=10, format="pdf")
fig.write_image(output_dir/f"enrichment_landscapes/network_{variable}.svg")
plot(fig, filename=str(output_dir/f"enrichment_landscapes/network_{variable}.html"), auto_open=False)
# %%
# Plot enrichment landscapes of microbiome phenotypes
for variable in list(oral_microbiome_genus.columns):
attribute = network_enrichment_scores[variable]
fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=True,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable=variable,
node_colormap="balance",
title=variable.split("genus_")[1],
titlefont_size=25,
annotation_text="",
color_range_min=-10,
color_range_max=10,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=False,
colorbar_annotation_text="",
height=500,
width=450,
)
pio.write_image(fig, output_dir/f"enrichment_landscapes_genera/network_{variable}.png", scale=10)
pio.write_image(fig, output_dir/f"enrichment_landscapes_genera/network_{variable}.pdf", scale=10, format="pdf")
fig.write_image(output_dir/f"enrichment_landscapes_genera/network_{variable}.svg")
plot(fig, filename=str(output_dir/f"enrichment_landscapes_genera/network_{variable}.html"), auto_open=False)
# %%
# Plot variable distribution on network of non-microbiome phenotypes
for variable in list(metadata.columns):
attribute = metadata_transformed[variable]
attribute_max = np.nanmax(attribute)
fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=False,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable=variable,
node_colormap="balance",
title=variable_styling_dict[variable].split(" ")[0],
titlefont_size=25,
annotation_text="",
color_range_min=attribute_max * -1,
color_range_max=attribute_max,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=True,
colorbar_annotation_text="",
height=500,
width=500,
)
pio.write_image(fig, output_dir/f"variable_landscapes/network_{variable}.png", scale=10)
pio.write_image(fig, output_dir/f"variable_landscapes/network_{variable}.pdf", scale=10, format="pdf")
fig.write_image(output_dir/f"variable_landscapes/network_{variable}.svg")
plot(fig, filename=str(output_dir/f"variable_landscapes/network_{variable}.html"), auto_open=False)
# %%
# Plot variable distribution on network of microbiome phenotypes
for variable in list(oral_microbiome_genus.columns):
attribute = zscore(oral_microbiome_genus_transformed[variable], nan_policy="omit")
attribute_max = np.nanmax(attribute)
fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=False,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable=variable,
node_colormap="balance",
title=variable.split("genus_")[1],
titlefont_size=25,
annotation_text="",
color_range_min=attribute_max * -1,
color_range_max=attribute_max,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=True,
colorbar_annotation_text="",
height=500,
width=500,
)
pio.write_image(fig, output_dir/f"variable_landscapes_genera/network_{variable}.png", scale=10)
pio.write_image(fig, output_dir/f"variable_landscapes_genera/network_{variable}.pdf", scale=10, format="pdf")
fig.write_image(output_dir/f"variable_landscapes_genera/network_{variable}.svg")
plot(fig, filename=str(output_dir/f"variable_landscapes_genera/network_{variable}.html"), auto_open=False)
# %%
# Plot colorbar to use for figures
cbar_fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=True,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable=variable,
node_colormap="balance",
title=variable_styling_dict[variable].split(" ")[0],
titlefont_size=25,
annotation_text="",
color_range_min=-10,
color_range_max=10,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=True,
colorbar_annotation_text="",
height=500,
width=500,
)
fig.write_image(output_dir/f"enrichment_landscapes/cbar_fig.svg")
# %% [markdown]
# # Cross-correlation matrices of network enrichment scores
# %%
# Microbiome phenotypes
import dash_bio
genus_variables = [col for col in network_enrichment_scores if "genus" in col]
network_enrichment_scores_corr = network_enrichment_scores[genus_variables].corr(method="spearman")
network_enrichment_scores_corr = network_enrichment_scores_corr
labels = [" ".join(var.split("_")[1:]) for var in list(network_enrichment_scores_corr.index)]
plot = dash_bio.Clustergram(
data=network_enrichment_scores_corr,
column_labels=labels,
row_labels=labels,
height=1600,
width=1700,
color_map="balance_r",
cluster="all",
center_values=False,
link_method="ward"
)
heatmap_trace=plot.data[-1]
heatmap_trace.update(colorbar_title='Spearman ρ ')
heatmap_trace.update(colorbar_xpad= 160)
plot_dict = plot.to_dict()
plot.write_html(output_dir/"enrichment_analysis/genera_cross_correlation_matrix.html")
plot.write_image(output_dir/"enrichment_analysis/genera_cross_correlation_matrix.png", format="png", scale=10)
plot.write_image(output_dir/"enrichment_analysis/genera_cross_correlation_matrix.svg", format="svg")
# %% [markdown]
# ## Metadata
# %%
variable_styling_dict_reverse = {v:k for k,v in variable_styling_dict.items()}
signif_envfit = compared_table[compared_table["envfit_p_fdr"] < 0.05].index.tolist()
signif_envfit = [variable_styling_dict_reverse[idx] for idx in signif_envfit]
corr_df = network_enrichment_scores.drop(oral_microbiome_genus, axis=1)
corr_df = corr_df.corr(method="spearman")
# %%
# Non-microbiome phenotypes
labels = [variable_styling_dict[var].split(" ")[0] for var in list(corr_df.index)]
plot = dash_bio.Clustergram(
data=corr_df,
column_labels=labels,
row_labels=labels,
height=800,
width=1100,
color_map="balance_r",
center_values=False,
)
# Adding heatmap to the clustergram layout
heatmap_trace=plot.data[-1]
heatmap_trace.update(colorbar_title='Spearman ρ ')
heatmap_trace.update(colorbar_xpad= 160)
plot.write_html(output_dir/"enrichment_analysis/metadata_cross_correlation_matrix.html")
plot.write_image(output_dir/"enrichment_analysis/metadata_cross_correlation_matrix.png", format="png", scale=10)
plot.write_image(output_dir/"enrichment_analysis/metadata_cross_correlation_matrix.svg", format="svg")
plot.show()
# %% [markdown]
# # Network clustering
# %%
# Perform KMeans clustering
from sklearn.cluster import KMeans
import numpy as np
positions = pd.DataFrame(nx.get_node_attributes(graph, "pos")).T
positions.columns = ["0", "1"]
clustering_input = positions.copy()
clustering_input.columns = [str(idx) for idx in list(range(clustering_input.shape[1]))]
n_clusters = 2
clustering = KMeans(n_clusters=2, random_state=42).fit(clustering_input)
positions["cluster"] = clustering.labels_
# %%
# Plot network with clustering annotation
attribute = positions["cluster"]
fig = plot_enrichment_landscape(
graph=graph,
attribute=attribute,
fade_nonsignificant_nodes=False,
network_enrichment_scores_signif=network_enrichment_scores_signif,
variable="demographics_age", # dummy variable to make the code work
node_colormap="RdBu_r",
title="",
annotation_text="",
color_range_min=0,
color_range_max=1,
nonsignif_opacity = 0.4,
node_line_color="black",
node_line_width=0.3,
show_colorbar=False,
colorbar_annotation_text="",
height=500,
width=450,
)
fig.write_html(output_dir/"cluster_analysis/graph_cluster.html")
fig.write_image(output_dir/"cluster_analysis/graph_cluster.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/graph_cluster.svg", format="svg")
fig.show()
# %% [markdown]
# # Perform subject-level group comparison
# %% [markdown]
# ### Non-microbiome phenotypes
# %%
# Transfer clustering information from node to subject level
import itertools
node_subject_mapping_idx_dict = {node:list(graph.nodes[idx]["sample"]) for idx,node in enumerate(graph.nodes)}
node_subject_mapping_dict = {node:list(graph.nodes[idx]["sample_names"]) for idx,node in enumerate(graph.nodes)}
all_subject_indices = sorted(set(itertools.chain(*node_subject_mapping_dict.values())))
node_subject_df = pd.DataFrame(0, index = all_subject_indices, columns = list(graph.nodes))
for node, subjects in node_subject_mapping_dict.items():
for subject in subjects:
node_subject_df.loc[subject, node] = 1
node_subject_df = node_subject_df.loc[metadata.index[metadata.index.isin(node_subject_df.index)]]
subject_group_df = node_subject_df.T.join(positions["cluster"]).groupby("cluster").sum().T
def determine_cluster(row):
if row[0] > 0 and row[1] > 0:
return -1
elif row[0] > 0:
return 0
elif row[1] > 0:
return 1
else:
return np.nan
subject_group_df["cluster"] = subject_group_df.apply(determine_cluster, axis=1)
subject_group_df.to_csv(output_dir/"cluster_analysis/subject_clustering.csv")
# %%
statistics_df = subject_group_df["cluster"].to_frame().join(metadata_plus_imaging).join(oral_microbiome_genus)
print("n of overlapping subjects:" , len(statistics_df[statistics_df["cluster"] == -1].index))
statistics_df = statistics_df[statistics_df["cluster"] != -1]
binary_variables = [col for col in statistics_df.columns if len(statistics_df[col].unique()) == 2]
continuous_variables = [col for col in statistics_df.columns if col not in binary_variables]
# %%
# Perform group comparison of non-microbiome phenotypes
import pingouin as pg
from scipy.stats import zscore
confound_variables = [
"demographics_age",
"demographics_sex",
"demographics_education_isced",
] + cardiovascular_risk_factors
def perform_linreg(results_df, dependent_variable):
df = statistics_df.copy()
df[confound_variables + [dependent_variable]] = zscore(df[confound_variables + [dependent_variable]], axis=0, nan_policy="omit")
model = pg.linear_regression(X=df[["cluster"] + confound_variables], y=df[dependent_variable], remove_na=True)
results_df.loc[dependent_variable,"p"] = model["pval"].values[1]
results_df.loc[dependent_variable,"coef"] = model["coef"].values[1]
results_df.loc[dependent_variable,"r2"] = model["r2"].values[1]
results_df.loc[dependent_variable,"CI[2.5%]"] = model["CI[2.5%]"].values[1]
results_df.loc[dependent_variable,"CI[97.5%]"] = model["CI[97.5%]"].values[1]
# %%
results_df = pd.DataFrame()
for dv in reversed(paro_variables + cognitive_scores+ neuropsychiatric_scores + imaging_means + inflammation + diet_scores): #
perform_linreg(results_df, dv)
results_df.index.rename("dv", inplace=True)
results_df["dv"] = results_df.index
results_df["p_fdr"] = pg.multicomp(results_df["p"].values, method="fdr_bh")[1]
# %%
gc_signif_variables = [
'paro_cal_mean',
'paro_dmft',
'paro_plaqueindex',
'paro_bop',
'cognition_g_factor_inverted',
'cognition_tmt_b_inverted',
'cognition_animal_naming_test',
'cognition_mini_mental_state_exam',
'imaging_thickness_volume_mean',
'inflammation_leukocytes',
]
# %%
import plotly.graph_objects as go
import pandas as pd
from scipy.stats import zscore
between = "cluster"
# Variables and Data Preparation
plotting_variables = gc_signif_variables
statistics_df_z = statistics_df[[between] + plotting_variables].copy()
statistics_df_z[plotting_variables] = statistics_df_z[plotting_variables].apply(lambda x: zscore(x, nan_policy="omit"), axis=0)
df = pd.melt(statistics_df_z, id_vars=["cluster"], value_vars=plotting_variables)
df["variable_styled"] = [variable_styling_dict[var] for var in df["variable"]]
# Reverse the order of variables
reversed_variables = list(reversed(df['variable_styled'].unique()))
# Create the Figure
fig = go.Figure()
# Add Box Plots with offsets in y to avoid overlapping
y_offset = 0.2
for i, var in enumerate(reversed_variables):
fig.add_trace(go.Box(x=df['value'][(df[between] == 1) & (df['variable_styled'] == var)],
y=[i - y_offset] * len(df[(df[between] == 1) & (df['variable_styled'] == var)]),
name='B',
marker_color='#6ca1bc',
orientation='h',
boxpoints=False,
showlegend=(i == 0))) # Show legend only for the first variable
fig.add_trace(go.Box(x=df['value'][(df[between] == 0) & (df['variable_styled'] == var)],
y=[i + y_offset] * len(df[(df[between] == 0) & (df['variable_styled'] == var)]),
name='A',
marker_color='#b35e4b', # Darker red color code
orientation='h',
boxpoints=False,
showlegend=(i == 0))) # Show legend only for the first variable
variable = df[df['variable_styled'] == var]['variable'].unique()[0]
coef = results_df.loc[variable, "coef"]
pval = results_df.loc[variable, "p_fdr"]
if pval < 0.001:
pval_styled = ' < 0.001'
else:
pval_styled = f' = {pval:.3f}'
fig.add_annotation(x=4.5, y=i,
text=f"βstd = {coef:.2f}", # p{pval_styled}
showarrow=False)
if pval < 0.001: asterisk = "***"
elif pval < 0.01: asterisk = "**"
elif pval < 0.05: asterisk = "*"
else: asterisk = ""
fig.add_annotation(x=2.5, y=i-0.05,
text=asterisk,
showarrow=False,
font=dict(size=18))
# Update Layout
fig.update_layout(
template="plotly_white",
title='',
xaxis_title="z",
yaxis_title="",
xaxis=dict(range=[-3, 4.5],
tickvals=[-2.5, 0, 2.5]), # Set x-axis limits
yaxis=dict(
tickvals=list(range(len(reversed_variables))),
ticktext=reversed_variables,
),
width=480,
height=900,
legend=dict(
title="",
x=-12, # Adjust the x position
y=-10, # Adjust the y position
xanchor="left", # Horizontal anchor at 'left'
yanchor="middle", # Vertical anchor at 'middle'
)
)
fig.write_html(output_dir/"cluster_analysis/group_comparison_brain_health_box.html")
fig.write_image(output_dir/"cluster_analysis/group_comparison_brain_health_box.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/group_comparison_brain_health_box.svg", format="svg")
fig.show()
# %% [markdown]
# ### Confound analysis
# %%
confound_variables = ['demographics_age',
'demographics_sex',
'demographics_education_isced',
'cvrisk_systolic_blood_pressure_mmhg',
'cvrisk_diastolic_blood_pressure_mmhg',
'cvrisk_BMI',
'cvrisk_smoking_currently',
'blood_cholesterol_mg_dl',
'blood_hdl_mg_dl',
'blood_ldl_mg_dl',
'blood_triglycerides_mg_dl',
'blood_hba1c']
# %%
def perform_linreg(results_df, dependent_variable, covariates):
df = statistics_df.copy()
if covariates == None:
df[dependent_variable] = zscore(df[dependent_variable], axis=0, nan_policy="omit")
model = pg.linear_regression(X=df["cluster"], y=df[dependent_variable], remove_na=True)
elif len(covariates) == 1:
df[covariates + [dependent_variable]] = zscore(df[covariates + [dependent_variable]], axis=0, nan_policy="omit")
model = pg.linear_regression(X=df[["cluster"] + covariates], y=df[dependent_variable], remove_na=True)
elif len(covariates) > 1:
df[covariates + [dependent_variable]] = zscore(df[covariates + [dependent_variable]], axis=0, nan_policy="omit")
model = pg.linear_regression(X=df[["cluster"] + covariates], y=df[dependent_variable], remove_na=True)
results_df.loc[dependent_variable,"p"] = model["pval"].values[1]
results_df.loc[dependent_variable,"coef"] = model["coef"].values[1]
results_df.loc[dependent_variable,"r2"] = model["r2"].values[1]
results_df.loc[dependent_variable,"CI[2.5%]"] = model["CI[2.5%]"].values[1]
results_df.loc[dependent_variable,"CI[97.5%]"] = model["CI[97.5%]"].values[1]
return results_df
# %%
p_df = pd.DataFrame()
coef_df = pd.DataFrame()
r2_df = pd.DataFrame()
for idx, dv in enumerate(paro_variables + inflammation + cognitive_scores+ neuropsychiatric_scores + imaging_means + diet_scores + confound_variables):
confound_iteration = []
confound_analysis_df_iter = pd.DataFrame()
dv_styled = variable_styling_dict[dv]
if " " in dv_styled: dv_styled = dv_styled.split(" ")[0]
confound_analysis_df_iter = perform_linreg(confound_analysis_df_iter, dv, covariates=None)
p_df.loc[dv_styled,"Unadjusted"] = confound_analysis_df_iter.iloc[0,0]
coef_df.loc[dv_styled,"Unadjusted"] = confound_analysis_df_iter.iloc[0,1]
r2_df.loc[dv_styled,"Unadjusted"] = confound_analysis_df_iter.iloc[0,2]
for confounder in confound_variables:
if dv == confounder: continue
confounder_styled = variable_styling_dict[confounder]
confound_iteration.append(confounder)
confound_analysis_df_iter = pd.DataFrame()
confound_analysis_df_iter = perform_linreg(confound_analysis_df_iter, dv, covariates=confound_iteration)
p_df.loc[dv_styled,"+ " + confounder_styled] = confound_analysis_df_iter.iloc[0,0]
coef_df.loc[dv_styled,"+ " + confounder_styled] = confound_analysis_df_iter.iloc[0,1]
r2_df.loc[dv_styled,"+ " + confounder_styled] = confound_analysis_df_iter.iloc[0,2]
# %%
import plotly.graph_objects as go
annotation_df = p_df.applymap(lambda p: '*' if p < 0.05 else '')
fig = go.Figure(data=go.Heatmap(
z=coef_df.values,
x=coef_df.columns,
y=coef_df.index,
text=annotation_df.values,
texttemplate="%{text}",
textfont={"size": 18, "color": "darkred"},
colorscale='RdBu_r',
zmid=0,
colorbar=dict(title='βstd')
))
fig.update_layout(
title='',
xaxis_title='',
yaxis_title='',
xaxis=dict(tickangle=-45),
yaxis=dict(autorange='reversed'),
width=900,
height=800,
margin=dict(l=250, r=50, b=150, t=50)
)
fig.show()
fig.write_html(output_dir/"cluster_analysis/confound_analysis.html")
fig.write_image(output_dir/"cluster_analysis/confound_analysis.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/confound_analysis.svg", format="svg")
# %% [markdown]
# ### Genera
# %%
from scipy.stats import zscore
from skbio.stats.composition import clr
microbiome_non_zero = oral_microbiome_genus + 1
clr_data = clr(microbiome_non_zero)
clr_df = pd.DataFrame(clr_data,
index=oral_microbiome_genus.index,
columns=oral_microbiome_genus.columns)
microbiome_z_columns = []
for dv in clr_df.columns:
z_col_name = f"{dv}_z"
statistics_df[z_col_name] = zscore(clr_df[dv], axis=0, nan_policy="omit")
microbiome_z_columns.append(z_col_name)
# %%
results_df = pd.DataFrame()
for dv in microbiome_z_columns:
perform_linreg(results_df, dv)
results_df.index.rename("dv", inplace=True)
results_df["dv"] = results_df.index
results_df["p_fdr"] = pg.multicomp(results_df["p"].values, method="fdr_bh")[1]
# %%
highest_change_pos = results_df[(results_df["p_fdr"]< 0.05) & (results_df["coef"]> 0)].sort_values(by="coef", ascending=False).iloc[:15].index.to_list()
highest_change_neg = results_df[(results_df["p_fdr"]< 0.05) & (results_df["coef"]< 0)].sort_values(by="coef").iloc[:15].index.to_list()
highest_change_pos = [var.split("_z")[0] for var in highest_change_pos]
highest_change_neg = [var.split("_z")[0] for var in highest_change_neg]
# %%
genera_signif = results_df[results_df["p_fdr"]< 0.05].sort_values(by="coef", ascending=False).index.to_list()
genera_signif = [var.split("_z")[0] for var in genera_signif]
# %%
import pandas as pd
import plotly.graph_objects as go
import plotly.express as px
# Filter the top 15 positive and top 15 negative coefficients
df_sorted = results_df.sort_values('coef', ascending=True)
df_top_pos = df_sorted.head(15)
df_top_neg = df_sorted.tail(15)
df_top = pd.concat([df_top_pos, df_top_neg])
colorscale = px.colors.qualitative.Prism + px.colors.qualitative.Plotly + px.colors.qualitative.Pastel # Choose your desired colorscale
# Define significance levels
def significance_indicator(p_value):
if p_value < 0.001:
return '***'
elif p_value < 0.01:
return '**'
elif p_value < 0.05:
return '*'
else:
return ''
df_top['significance'] = df_top['p'].apply(significance_indicator)
df_top['genus'] = df_top['dv'].apply(lambda x: " ".join(x.split('_')[1:]).split(" z")[0]) # Simplify genus names if necessary
# Create the bar plot
fig = go.Figure()
# Add bars
fig.add_trace(go.Bar(
x=df_top['genus'],
y=df_top['coef'],
error_y=dict(
type='data',
symmetric=False,
array=df_top['CI[97.5%]'] - df_top['coef'],
arrayminus=df_top['coef'] - df_top['CI[2.5%]']
),
marker=dict(color=df_top['coef'], colorscale='RdBu'), # dict(color=colorscale),#
hoverinfo='x+y',
))
# Add significance annotations
for i, row in df_top.iterrows():
if row['significance']:
fig.add_annotation(y=row["coef"] + 0.4 if row['coef'] > 0 else row["coef"] -0.4, x=row['genus'],
text=row['significance'],
showarrow=False, font=dict(color='black'), textangle=0)
fig.add_annotation(x=2.5, y=1.5,
text="*:pFDR<0.05 **:pFDR<0.01 ***:pFDR<0.001",
showarrow=False, font=dict(color='black'))
# Horizontal line at y=0
fig.add_hline(y=0, line=dict(color="grey", width=1))
# Update layout
fig.update_layout(
title="",
xaxis_title="",
yaxis_title="βstd",
template="plotly_white",
height=400,
width=850,
xaxis=dict(tickangle=45) # Rotate the x-ticks by -45 degrees
)
fig.write_html(output_dir/"cluster_analysis/abundancy_cluster_barplot_coefficients.html")
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_barplot_coefficients.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_barplot_coefficients.svg", format="svg")
# Show the plot
fig.show()
# %%
df_sorted = results_df.sort_values('coef', ascending=True)
df_top = df_sorted
df_top['genus'] = df_top['dv'].apply(lambda x: " ".join(x.split('_')[1:]).split(" z")[0]) # Simplify genus names if necessary
# Create the horizontal bar plot
fig = go.Figure()
# Add bars with orientation set to horizontal
fig.add_trace(go.Bar(
y=df_top['genus'], # Use y for horizontal bar plot
x=df_top['coef'], # Use x for horizontal bar plot
error_x=dict( # Adjust error_x instead of error_y for horizontal errors
type='data',
symmetric=False,
array=df_top['CI[97.5%]'] - df_top['coef'],
arrayminus=df_top['coef'] - df_top['CI[2.5%]'],
width=1
),
marker=dict(color=df_top['coef'], colorscale='RdBu'), # Update marker color
hoverinfo='y+x', # Adjust hover info for horizontal orientation
orientation='h' # Set orientation to horizontal
))
# Vertical line at x=0
fig.add_vline(x=0, line=dict(color="grey", width=1))
# Update layout for horizontal bar plot
fig.update_layout(
title="",
yaxis_title="", # Update y-axis title for horizontal display
xaxis_title="βstd",
template="plotly_white",
height=1400, # Adjust height for horizontal plot
width=700, # Adjust width for horizontal plot
)
fig.write_html(output_dir/"cluster_analysis/abundancy_cluster_barplot_coefficients_all.html")
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_barplot_coefficients_all.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_barplot_coefficients_all.svg", format="svg")
# Save or display the horizontal bar plot
fig.show()
# %%
import plotly.graph_objects as go
import pandas as pd
import plotly.colors as pcolors # Import Plotly color scales
import plotly.express as px
# Define a custom color palette (e.g., using Plotly's predefined color scales)
colorscale = px.colors.qualitative.Prism + px.colors.qualitative.Plotly + px.colors.qualitative.Pastel # Choose your desired colorscale
genera_cols = genera_signif
# Copying and transforming the dataset
plotting_df = statistics_df[["cluster"] + highest_change_pos + highest_change_neg].copy()
genera_cols_styled = ["_".join(idx.split("_")[1:]) for idx in highest_change_pos + highest_change_neg]
plotting_df.columns = ["cluster"] + genera_cols_styled
plotting_df[genera_cols_styled] = plotting_df[genera_cols_styled].div(statistics_df[oral_microbiome_genus.columns].sum(axis=1), axis=0)
# Melting the dataframe for long-form plotting
plotting_df_0 = plotting_df[plotting_df["cluster"] == 0].mean(axis=0)
plotting_df_1 = plotting_df[plotting_df["cluster"] == 1].mean(axis=0)
# Sort each cluster's values in descending order
plotting_df_0 = plotting_df_0.sort_values(ascending=False)
plotting_df_1 = plotting_df_1.sort_values(ascending=False)
# Combine into a DataFrame
plotting_df = pd.concat([plotting_df_0, plotting_df_1], axis=1).T
plotting_df = plotting_df * 100
plotting_df["cluster"] = [0, 1]
# Reorder columns based on sorted indexes for cluster 0 and cluster 1
sorted_columns = plotting_df_0.drop("cluster").index.tolist()
# Creating the main Figure
fig = go.Figure()
# Looping through each genus in sorted order to add a bar trace
color_idx = 0
for genus in sorted_columns:
color = colorscale[color_idx % len(colorscale)] # Get color from colorscale
color_idx += 1
fig.add_trace(go.Bar(
x=plotting_df["cluster"],
y=plotting_df[genus],
name=genus,
text=genus if plotting_df[genus].max() > 1 else "", # Only show text if the max value is greater than 1
hovertemplate=f"% {genus} " + "%{y:.3f}<extra></extra>", # Display genus name on hover
marker_color=color, # Assign color to the trace
))
# Updating layout to match the original requirements
fig.update_layout(
template="plotly_white",
height=600,
width=550,
barmode='stack',
xaxis_title="Cluster",
yaxis_title="Abundance (%)",
title="",
showlegend=True,
)
# Customizing x-axis to show only 0 and 1
fig.update_xaxes(
tickvals=[0, 1],
ticktext=["A", "B"]
)
# save html
fig.write_html(output_dir/"cluster_analysis/abundancy_cluster_high_diff_genera.html")
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_high_diff_genera.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_high_diff_genera.svg", format="svg")
# Showing the plot
fig.show()
# %%
import plotly.graph_objects as go
import pandas as pd
import plotly.colors as pcolors # Import Plotly color scales
import plotly.express as px
# Define a custom color palette (e.g., using Plotly's predefined color scales)
colorscale = px.colors.qualitative.Prism + px.colors.qualitative.Plotly + px.colors.qualitative.Pastel # Choose your desired colorscale
genera_cols = genera_signif
# Copying and transforming the dataset
plotting_df = statistics_df[["cluster"] + genera_cols].copy()
genera_cols_styled = ["_".join(idx.split("_")[1:]) for idx in genera_cols]
plotting_df.columns = ["cluster"] + genera_cols_styled
plotting_df[genera_cols_styled] = plotting_df[genera_cols_styled].div(plotting_df[genera_cols_styled].sum(axis=1), axis=0)
# Melting the dataframe for long-form plotting
plotting_df_0 = plotting_df[plotting_df["cluster"] == 0].mean(axis=0)
plotting_df_1 = plotting_df[plotting_df["cluster"] == 1].mean(axis=0)
# Sort each cluster's values in descending order
plotting_df_0 = plotting_df_0.sort_values(ascending=False)
plotting_df_1 = plotting_df_1.sort_values(ascending=False)
# Combine into a DataFrame
plotting_df = pd.concat([plotting_df_0, plotting_df_1], axis=1).T
plotting_df = plotting_df * 100
plotting_df["cluster"] = [0, 1]
# Reorder columns based on sorted indexes for cluster 0 and cluster 1
sorted_columns = plotting_df_0.drop("cluster").index.tolist()
# Creating the main Figure
fig = go.Figure()
# Looping through each genus in sorted order to add a bar trace
color_idx = 0
for genus in sorted_columns:
color = colorscale[color_idx % len(colorscale)] # Get color from colorscale
color_idx += 1
fig.add_trace(go.Bar(
x=plotting_df["cluster"],
y=plotting_df[genus],
name=genus,
text=genus if plotting_df[genus].max() > 2 else "", # Only show text if the max value is greater than 1
hovertemplate=f"% {genus} " + "%{y:.3f}<extra></extra>", # Display genus name on hover
marker_color=color, # Assign color to the trace
))
# Updating layout to match the original requirements
fig.update_layout(
template="plotly_white",
height=900,
width=600,
barmode='stack',
xaxis_title="Cluster",
yaxis_title="Abundance (%)",
title="",
showlegend=False,
)
# Customizing x-axis to show only 0 and 1
fig.update_xaxes(
tickvals=[0, 1],
ticktext=["A", "B"]
)
# save html
fig.write_html(output_dir/"cluster_analysis/abundancy_cluster_all_genera.html")
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_all_genera.png", format="png", scale=10)
fig.write_image(output_dir/"cluster_analysis/abundancy_cluster_all_genera.svg", format="svg")
# Showing the plot
fig.show()
# %% [markdown]
# # Dominant genera
# %%
# modified from tmap documentation https://tmap.readthedocs.io/en/latest/
import plotly.graph_objects as go
from collections import Counter
import numpy as np
import plotly.subplots
import plotly.express as px
def plot_dominance_landscape(
graph=None,
dominance_df=None,
thresh_top_enrichment=5, # min amount of top enrichment nodes a variable must have to plot it
colormap=px.colors.qualitative.Prism, # Using a Plotly colorscale
node_line_color="black",
node_line_width=0.3, # node line width
opacity=0.9, # node opacity
width=900,
height=700,
):
# Assuming 'pos' is an attribute in the graph nodes that stores position
node_pos = np.array([graph.nodes[node]['pos'] for node in graph.nodes()])
# Edges for the graph
xs, ys = [], []
for edge in graph.edges:
xs += [node_pos[edge[0]][0], node_pos[edge[1]][0], None]
ys += [node_pos[edge[0]][1], node_pos[edge[1]][1], None]
fig = plotly.subplots.make_subplots(rows=1, cols=1)
# Add edges to the plot
fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines",
line=dict(width=1, color="#8E9DA2"), showlegend=False), row=1, col=1)
# Calculate which feature is most dominant per node
feature_indices = np.argmax(dominance_df.values, axis=1)
tmp = [dominance_df.columns[index] for index in feature_indices]
# Count the number of top enrichment nodes per variable
t = Counter(tmp)
enrichment_features = {fea for fea, count in t.items() if count >= thresh_top_enrichment}
cmap = {feature: colormap[i % len(colormap)] for i, feature in enumerate(enrichment_features)}
# Plot nodes based on the dominant feature and use colormap for nodes
for feature in enrichment_features:
indices = [i for i, f in enumerate(tmp) if f == feature]
fig.add_trace(go.Scatter(
x=node_pos[indices, 0], y=node_pos[indices, 1], mode='markers',
marker=dict(size=15, color=cmap[feature], opacity=opacity, line=dict(width=node_line_width, color=node_line_color)),
name=f'{variable_styling_dict[feature].split(" ")[0] if feature in variable_styling_dict.keys() else feature.replace("_"," ")} ({t[feature]})',
showlegend=True
), row=1, col=1)
# Update layout
fig.update_layout(
width=width, height=height, hovermode='closest', plot_bgcolor="white", paper_bgcolor="white",
xaxis=dict(showgrid=False, zeroline=False, showticklabels=False),
yaxis=dict(showgrid=False, zeroline=False, showticklabels=False)
)
return fig
# %%
dominance_df = network_enrichment_scores[oral_microbiome_genus.columns]
colormap = plotly.colors.qualitative.Prism + plotly.colors.qualitative.Plotly
dominance_df.columns = [" ".join(col.split("_")[1:]) for col in dominance_df.columns]
fig = plot_dominance_landscape(
graph=graph,
dominance_df=dominance_df,
colormap=colormap,
thresh_top_enrichment=5,
node_line_color="black",
node_line_width=0.3,
opacity=0.9,
width=800,
height=700,
)
fig.write_html(output_dir/"enrichment_analysis/dominance_all_genera.html")
fig.write_image(output_dir/"enrichment_analysis/dominance_all_genera.svg")
fig.write_image(output_dir/"enrichment_analysis/dominance_all_genera.png", format="png", scale=10)
fig.show()

analysis.ipynb at commit 5fee720, no license · at the source

Overview

Authors: Marvin Petersen¹, Carolin Walther², Katrin Borof², Guido Heydecke³, Thomas Beikler², Malik Alawi⁴, Christian Müller⁴, Felix L. Nägele¹, Birgit-Christiane Zyriax⁵, Jens Fiehler⁶, Jürgen Gallinat⁷, Simone Kühn^7,8, Raphael Twerenbold^9,10,11,12, Corinna Bang¹³, Götz Thomalla¹, Bastian Cheng¹, Ghazal Aarabi²

ORCID iDs: Marvin Petersen, Bastian Cheng, Ghazal Aarabi

13 affiliations

Department of Neurology, University Medical Centre Hamburg-Eppendorf, Hamburg, Germany
Department of Periodontics, Preventive and Restorative Dentistry, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
Department of Prosthetic Dentistry, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
Bioinformatics Core, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
Midwifery Science-Health Services Research and Prevention, Institute for Health Services Research in Dermatology and Nursing (IVDP), University Medical Center Hamburg-Eppendorf, Hamburg, Germany
Department of Neuroradiology, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
Department of Psychiatry and Psychotherapy, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
Lise Meitner Group for Environmental Neuroscience, Max Planck Institute for Human Development, Berlin, Germany
Department of Cardiology, University Heart and Vascular Center, Hamburg, Germany
Epidemiological Study Center, University Medical Center Hamburg-Eppendorf, Hamburg, Germany
German Center for Cardiovascular Research (DZHK), Partner Site Hamburg/Kiel/Luebeck, Hamburg, Germany
University Center of Cardiovascular Science, University Heart and Vascular Center, Hamburg, Germany
Institute of Clinical Molecular Biology, Kiel University, Kiel, Germany

Institutions: University Medical Center Hamburg-Eppendorf (Germany); Max Planck Institute for Human Development (Germany); German Centre for Cardiovascular Research (Germany); Christian-Albrechts-Universität zu Kiel (Germany)

Journal: EBioMedicine, volume 128, article 106312

Dates: received 28 December 2025; accepted 12 May 2026; published online 30 May 2026; in print June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1016/j.ebiom.2026.106312 · PMID 42217285 · PMCID PMC13241667 · OpenAlex W7162878770

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), cognitive (subfield)

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

Keywords: Oral microbiome, Periodontitis, Oral-brain axis, Brain health, Cognition, Neuroimaging

MeSH: Brain*, Gingiva*, Microbiota*, Periodontitis*, Aged, Bacteria, Biomarkers, Cognition, Female, Humans, Male, Middle Aged, RNA, Ribosomal, 16S (* major topic)

Topic: Gut microbiota and health (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: University Medical Center Hamburg-Eppendorf

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

Abstract

Background: Periodontitis has gained attention as a key factor associated with cognitive decline and Alzheimer's dementia. However, the relationship between periodontitis-related oral microbiota shifts and brain health in the general population remains unclear.

Methods: We investigated the oral microbiome–brain axis in 1026 participants from the population-based PAROMIND Study. Using 16S rRNA gene amplicon sequencing of subgingival crevicular fluid, we inferred via topological data analysis a microbiota similarity network. This network, which distills the complex high-dimensional data into an interpretable map of microbial similarity, revealed a continuous disease gradient mirroring the microbial pathogenicity spectrum, from taxa of low periodontal pathogenicity (e.g., Streptococcus) to periodontitis-associated taxa (e.g., Porphyromonas, Fusobacterium). Leveraging this network, we systematically examined associations between periodontal microbiota profiles and 40 brain health-related phenotypes, including cognition, brain structure, mental health, inflammatory biomarkers, diet, vascular risk factors, and demographics.

Findings: Higher abundance of periodontitis-related bacterial taxa was associated with poorer cognitive performance, elevated leucocyte counts, and lower MIND diet adherence after covariate adjustment. Complementary forward model selection analysis supported the links to cognitive performance and inflammation, and additionally identified a significant association with brain structure (cortical thickness and subcortical volume). We identified associations with both established genera (Porphyromonas) and taxa not previously implicated in brain health (Fretibacterium, Tannerella, Dialister).

Interpretation: These findings from a large cohort advance the understanding of the oral microbiome–brain axis, highlighting specific microbial profiles linked to subclinical cognitive, structural, and inflammatory brain health markers. By demonstrating these links in a non-demented population, our study suggests that monitoring the oral microbiome could inform early risk assessment for cognitive decline, positioning periodontal health as an accessible target for early intervention strategies.

Funding: Deutsche Forschungsgemeinschaft.

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

Repository

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

csi-hamburg/oral_microbiome_brain_health

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 5fee7207cc6aa54a16c588f93bb9b6f9ad418497, 1 June 2026

Languages: Jupyter (1)

Size: 2 files, 1 script

Software Heritage: not archived

Found in: “Data sharing statement”

Holds: README, 1 notebook

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

Tools: NetworkX (1 file), NumPy (1 file), pandas (1 file), Pingouin (1 file), Plotly (1 file), rpy2 (1 file), scikit-learn (1 file), SciPy (1 file), statsmodels (1 file)

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

28 September 2026: the link answers

2 files

analysis.ipynb, Jupyter, 2,057 lines, 4 matches
README.md, Text, 2 lines

At the source: github.com/csi-hamburg/oral_microbiome_brain_health

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:

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

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

Data

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

Data sharing statement

Sequencing data generated during this study have been deposited as FASTQ files in the European Nucleotide Archive (ENA) under accession number PRJEB89258. Corresponding non-microbiome phenotype data from the HCHS are not publicly available due to data protection policies that ensure participant confidentiality. However, this data is available to qualified researchers upon reasonable request to the HCHS steering committee. The analysis code for this work is publicly available on GitHub (https://github.com/csi-hamburg/oral_microbiome_brain_health). Interactive versions of the plots can be found on OSF (https://osf.io/vqj8m/).

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

Versions

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

Version 2, 28 September 2026

Authors: added Bastian Cheng (0000-0003-2434-1822); Ghazal Aarabi (0000-0001-5484-2594); removed Bastian Cheng; Ghazal Aarabi

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 17 authors, 6 keywords, 13 MeSH terms, 1 funder, 66 references.

Cite

This paper

Petersen, M., Walther, C., Borof, K., Heydecke, G., Beikler, T., Alawi, M., Müller, C., Nägele, F. L., Zyriax, B.-C., Fiehler, J., Gallinat, J., Kühn, S., Twerenbold, R., Bang, C., Thomalla, G., Cheng, B., & Aarabi, G. (2026). Subgingival microbiota composition is associated with brain health in the general population-the PAROMIND study. EBioMedicine, 128, 106312. https://doi.org/10.1016/j.ebiom.2026.106312

BibTeX

@article{petersen2026subgingival, author = {Petersen, Marvin and Walther, Carolin and Borof, Katrin and Heydecke, Guido and Beikler, Thomas and Alawi, Malik and Müller, Christian and Nägele, Felix L. and Zyriax, Birgit-Christiane and Fiehler, Jens and Gallinat, Jürgen and Kühn, Simone and Twerenbold, Raphael and Bang, Corinna and Thomalla, Götz and Cheng, Bastian and Aarabi, Ghazal}, title = {{Subgingival microbiota composition is associated with brain health in the general population-the PAROMIND study}}, journal = {EBioMedicine}, year = {2026}, month = may, volume = {128}, pages = {106312}, publisher = {Elsevier}, issn = {2352-3964}, doi = {10.1016/j.ebiom.2026.106312}, url = {https://doi.org/10.1016/j.ebiom.2026.106312}, pmid = {42217285}, pmcid = {PMC13241667} }

RIS

TY - JOUR AU - Petersen, Marvin AU - Walther, Carolin AU - Borof, Katrin AU - Heydecke, Guido AU - Beikler, Thomas AU - Alawi, Malik AU - Müller, Christian AU - Nägele, Felix L. AU - Zyriax, Birgit-Christiane AU - Fiehler, Jens AU - Gallinat, Jürgen AU - Kühn, Simone AU - Twerenbold, Raphael AU - Bang, Corinna AU - Thomalla, Götz AU - Cheng, Bastian AU - Aarabi, Ghazal TI - Subgingival microbiota composition is associated with brain health in the general population-the PAROMIND study T2 - EBioMedicine J2 - EBioMedicine PY - 2026 DA - 2026/05/30 VL - 128 SP - 106312 SN - 2352-3964 PB - Elsevier DO - 10.1016/j.ebiom.2026.106312 UR - https://doi.org/10.1016/j.ebiom.2026.106312 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.ebiom.2026.106312", "type": "article-journal", "title": "Subgingival microbiota composition is associated with brain health in the general population-the PAROMIND study", "container-title": "EBioMedicine", "author": [ { "family": "Petersen", "given": "Marvin" }, { "family": "Walther", "given": "Carolin" }, { "family": "Borof", "given": "Katrin" }, { "family": "Heydecke", "given": "Guido" }, { "family": "Beikler", "given": "Thomas" }, { "family": "Alawi", "given": "Malik" }, { "family": "Müller", "given": "Christian" }, { "family": "Nägele", "given": "Felix L." }, { "family": "Zyriax", "given": "Birgit-Christiane" }, { "family": "Fiehler", "given": "Jens" }, { "family": "Gallinat", "given": "Jürgen" }, { "family": "Kühn", "given": "Simone" }, { "family": "Twerenbold", "given": "Raphael" }, { "family": "Bang", "given": "Corinna" }, { "family": "Thomalla", "given": "Götz" }, { "family": "Cheng", "given": "Bastian" }, { "family": "Aarabi", "given": "Ghazal" } ], "container-title-short": "EBioMedicine", "volume": "128", "page": "106312", "DOI": "10.1016/j.ebiom.2026.106312", "PMID": "42217285", "PMCID": "PMC13241667", "ISSN": "2352-3964", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.ebiom.2026.106312", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 30 ] ] } }

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

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.