Quantifying extra-lesional interhemispheric cortical asymmetry in focal cortical dysplasia type II.

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 ↔ Bonn_visualization.ipynb, lines 449–497 · score 0.62 · occipital lobe, temporal lobe, frontal lobe, insular, subset, lesions
[2] § METHODOLOGY › Surface‐based morphometry ↔ Bonn_visualization.ipynb, lines 518–524 · score 0.55 · mri_vol2surf, lesion volume, cortical surface, space

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 · 5,297 lines · 195 KB · MIT · 2 matches

# %% [markdown]
# # Visualization of Dataset
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Define the age intervals and corresponding labels
age_bins = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
age_labels = [
"0-5 years", "6-10 years", "11-15 years", "16-20 years", "21-25 years",
"26-30 years", "31-35 years", "36-40 years", "41-45 years", "46-50 years",
"51-55 years", "56-60 years", "61-65 years"
]
# Drop missing values in the 'age_epilepsyonset' column
df2 = df2.dropna(subset=['age_epilepsyonset'])
# Bin the 'age_epilepsyonset' data into the defined age intervals
df2['age_interval'] = pd.cut(df2['age_epilepsyonset'], bins=age_bins, labels=age_labels, right=False)
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='age_interval', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Distribution of Epilepsy Onset Age", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("Age Interval", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(["Epilepsy Onset Age"], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Define the age intervals and corresponding labels
age_bins = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
age_labels = [
"0-5 years", "6-10 years", "11-15 years", "16-20 years", "21-25 years",
"26-30 years", "31-35 years", "36-40 years", "41-45 years", "46-50 years",
"51-55 years", "56-60 years", "61-65 years"
]
# Drop missing values in the 'age_scan' column
df2 = df2.dropna(subset=['age_scan'])
# Bin the 'age_scan' data into the defined age intervals
df2['age_interval'] = pd.cut(df2['age_scan'], bins=age_bins, labels=age_labels, right=False)
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='age_interval', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Distribution of age at scan onset", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("Age Interval", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(["age_scan"], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Define the age intervals and corresponding labels
age_bins = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
age_labels = [
"0-5 years", "6-10 years", "11-15 years", "16-20 years", "21-25 years",
"26-30 years", "31-35 years", "36-40 years", "41-45 years", "46-50 years",
"51-55 years", "56-60 years", "61-65 years"
]
# Drop missing values in the 'age_op' column
df2 = df2.dropna(subset=['age_op'])
# Bin the 'age_op' data into the defined age intervals
df2['age_interval'] = pd.cut(df2['age_op'], bins=age_bins, labels=age_labels, right=False)
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='age_interval', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Distribution of age at operation onset", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("Age Interval", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(["Age at operation onset"], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Define the age intervals and corresponding labels
age_bins = [1, 2, 3, 4, 5, 6, 7, 8]
age_labels = [
"less than 1 year", "1-2 years", "3-4 years", "5-6 years", "7-8 years",
"9-10 years", "11-12 years"]
# Drop missing values in the 'time_follow-up' column
df2 = df2.dropna(subset=['time_follow-up'])
# Bin the 'time_follow-up' data into the defined age intervals
df2['age_interval'] = pd.cut(df2['time_follow-up'], bins=age_bins, labels=age_labels, right=False)
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='age_interval', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Distribution of time_follow-up", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("Age Interval", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(["time_follow-up"], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Drop missing values in the 'group' column
df2 = df2.dropna(subset=['group'])
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='group', hue='group', shrink=0.8, multiple="dodge",
palette={"fcd": "teal", "hc": "cyan"}, # Define colors for each group
edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Number of Subjects with FCD and Healthy Controls", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("Group", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(['fcd','hc'], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Drop missing values in the 'histopathology' column
df2 = df2.dropna(subset=['histopathology'])
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='histopathology', hue='histopathology', shrink=0.8, multiple="dodge",
palette={"IIa": "teal", "IIb": "cyan"}, # Define colors for each group
edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Distrubution of patients with fcd", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("histopathology", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(['fcdIIa','fcdIIb'], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Assuming cdf_fcd contains your data
cdf_fcd = pd.DataFrame(df2)
# Filter the DataFrame for the 'fcd' group
fcd_data = cdf_fcd[cdf_fcd['group'] == 'fcd'].copy() # Use .copy() to avoid SettingWithCopyWarning
# Replace NaN with a string to handle them in counts
fcd_data.loc[:, 'histopathology'] = fcd_data['histopathology'].fillna('No operation')
# Replace 'NaN' with 'No Operation' in the 'histopathology' column
cdf_fcd['histopathology'] = cdf_fcd['histopathology'].replace('NaN', 'No Operation')
# If we have actual NaN values (as missing data) and want to replace them:
cdf_fcd['histopathology'] = cdf_fcd['histopathology'].fillna('No Operation')
# Count the number of patients in each histopathology category within the 'fcd' group
group_counts_fcd = fcd_data['histopathology'].value_counts()
# Define colors for categories if known
colors = ['teal', 'cyan', 'darkgreen'] # Adjust colors as needed for each category
# Plot the histogram for all histopathology categories in the 'fcd' group
group_counts_fcd.plot(kind='bar', color=colors)
# Adding titles and labels
plt.title('Distribution of Histopathology Types diagnosis in FCD Group', fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel('Histopathology Type', fontsize=14, fontweight='bold')
plt.ylabel('Number of Patients', fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Drop missing values in the 'sex' column
df2 = df2.dropna(subset=['sex'])
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='sex', hue='sex', shrink=0.8, multiple="dodge",
palette={"F": "teal", "M": "cyan"}, # Define colors for each sex
edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Distribution of sex", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("sex", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(['Female','Male'], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Define the angel levels and corresponding labels
angel_bins = [1, 2, 3, 4, 5, 6, 7]
angel_labels = [
"IA,Completely seizure free since surgery", "IB,Nondisabling simple partial seizures only since surgery", "IIA,Initially free of disabling seizures but has rare seizures now", "IIB,Rare disabling seizures since surgery", "IIIA,Engel IIIA: Worthwhile seizure reduction",
"IVB,No appreciable change"]
# Drop missing values in the '1year_outcome' column
df2 = df2.dropna(subset=['1year_outcome'])
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='1year_outcome', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("distribution of 1year outcome", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("1year outcome", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
# Load the full path
file_path = '/Volumes/groups/tohkagroup/Bonn_Epilepsy/derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv'
# === Check if the file exists before reading ===
if not os.path.exists(file_path):
raise FileNotFoundError(f"File not found: {file_path}")
# === Load the TSV file ===
df2 = pd.read_csv(file_path, sep='\t')
# === Drop missing values in the 'latest_outcome' column ===
df2 = df2.dropna(subset=['latest_outcome'])
# === Plotting the histogram ===
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='latest_outcome', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# === Customize the plot ===
plt.title("Distribution of Latest Outcome", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("Latest Outcome (Engel Classification)", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# === Display the plot ===
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Drop missing values in the 'hemisphere' column
df2 = df2.dropna(subset=['hemisphere'])
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='hemisphere', hue='hemisphere', shrink=0.8, multiple="dodge",
palette={"L": "teal", "R": "cyan"}, # Define colors for each hemisphere
edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("Hemisphere affected by the FCD", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("hemisphere", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Adding a legend
plt.legend(['Right','Left'], fontsize=12, loc="upper right", title="Legend", title_fontsize='13')
# Display the plot
plt.tight_layout()
plt.show()
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_csv('derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv', sep='\t')
# Define the lobe levels and corresponding labels
lobe_bins = [1, 2, 3, 4, 5, 6]
lobe_labels = [
"frontal lobe", "temporal lobe", "pariatal lobe", "occipital lobe", "insular lobe"]
# Drop missing values in the 'lobe' column
df2 = df2.dropna(subset=['lobe'])
# Plotting the histogram
plt.figure(figsize=(12, 7))
sns.histplot(data=df2, x='lobe', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Customizing the plot
plt.title("distribution of lobe affected by FCD", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("lobe levels", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Display the plot
plt.tight_layout()
plt.show()
# %% [markdown]
# # Distribution of affected lobe by lesion in Stage 3
# %%
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Load the participants.aparc.tsv data into a DataFrame
df2 = pd.read_excel('/Volumes/groups/tohkagroup/Bonn_Epilepsy/derivatives_fsaverage/freesurfer7.4.1/participants.aparc_stage3.xlsx')
# Define the lobe levels and corresponding labels
lobe_bins = [1, 2, 3, 4, 5, 6]
lobe_labels = [
"frontal lobe", "temporal lobe", "pariatal lobe", "occipital lobe", "insular lobe"]
# Drop missing values in the 'lobe' column
df2 = df2.dropna(subset=['lobe'])
# Total number of subjects for percentage calculation
total= len(df2)
# Plotting the histogram
plt.figure(figsize=(12, 7))
ax=sns.histplot(data=df2, x='lobe', shrink=0.8, color="teal", edgecolor="black", linewidth=1.2)
# Add percentage labels on top of each bar
for p in ax.patches:
count = p.get_height()
percentage = (count / total) * 100
ax.annotate(
f"{percentage:.1f}%",
(p.get_x() + p.get_width() / 2., count),
ha='center',
va='bottom',
fontsize=12,
fontweight='bold'
)
# Customizing the plot
plt.title("Distribution of lobe affected by FCD II", fontsize=18, fontweight='bold', color='darkblue')
plt.xlabel("lobe", fontsize=14, fontweight='bold')
plt.ylabel("Number of Subjects", fontsize=14, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle='--', alpha=0.7)
# Display the plot
plt.tight_layout()
plt.show()
# %% [markdown]
# ## Asymmetry Analysis
# %% [markdown]
# # Set the lesion mask to the same directory ( copy all the lesion mask to their own subject files)
# %% [markdown]
# # In high-performance servers like Kudos, perform the following three steps
# %% [markdown]
# Step 1: Convert NIfTI to Freesurfer Volume
# %%
mri_vol2vol --mov sub-00146_acq-T2sel_FLAIR_roi.nii.gz \
--targ /research/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/mri/orig.mgz \
--regheader \
--o sub-00146_lesion_fs.mgz
# %% [markdown]
# Step 2: Convert Lesion Volume to a Surface Label
# Once the lesion is in Freesurfer’s space, map it onto the cortical surface: (but please consider it based on your lesion lateralization use either right or left hemisphere ) for example for this subject the lesion is located in left hemisphere, and of course, export the directory to put the result in the actual directory:
# %%
export SUBJECTS_DIR=/research/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna
mri_vol2surf --mov sub-00146_lesion_fs.mgz --o sub-00146_lesion_surf.mgz --regheader sub-00146 --projfrac 0.5 --hemi rh
# %% [markdown]
# Step 3: Identify ROIs Overlapping with the Lesion
# Better Solution: Using mri_segstats to Get ROI Overlaps at Once
# %%
export SUBJECTS_DIR=/research/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna
mri_segstats --annot sub-00146 rh aparc \
--i $SUBJECTS_DIR/sub-00146/sub-00146_lesion_surf.mgz \
--sum $SUBJECTS_DIR/sub-00146/stats/rh.lesion_roi_overlap.txt
# %% [markdown]
# # Removing cortical overlap:
# %% [markdown]
# # Filtering overlapped ROIs from lh.w-g.pct.stats and rh.w-g.pct.stats Based on StructName
# %%
import os
import pandas as pd
from io import StringIO
# Input and output paths
input_file = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats/rh.lesion_roi_overlap.txt"
output_file = input_file.replace(".txt", "_filtered.txt")
# Read the original file
with open(input_file, "r") as f:
lines = f.readlines()
# Locate the '# ColHeaders' line
col_header_idx = None
for i, line in enumerate(lines):
if line.startswith("# ColHeaders"):
col_header_idx = i
break
if col_header_idx is None:
raise ValueError("Could not find '# ColHeaders' line in the file.")
# Extract column names from that line
column_names = lines[col_header_idx].strip().replace("# ColHeaders", "").split()
# Header + table sections
header_lines = lines[:col_header_idx + 1]
table_lines = lines[col_header_idx + 1:]
# Load table into DataFrame
df = pd.read_csv(StringIO("".join(table_lines)), delim_whitespace=True, names=column_names)
# Filter for non-zero Mean
filtered_df = df[df["Mean"] != 0]
# Save back the filtered data
with open(output_file, "w") as out:
# Write header first
for line in header_lines:
out.write(line)
# Write filtered table without header row (ColHeaders already written)
filtered_df.to_csv(out, sep="\t", index=False, header=False)
print(f"Filtered data saved to: {output_file}")
# %% [markdown]
# # Removing filtered ROIs from stats file for all measurements
# %%
import pandas as pd
from io import StringIO
#Input files
lh_stats_file = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats/lh.w-g.pct.stats"
rh_overlap_file = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats/rh.lesion_roi_overlap_filtered.txt"
output_file = lh_stats_file.replace(".stats", ".filtered.stats")
#Step 1: Read the overlap file and collect Index values
indices_to_exclude = set()
with open(rh_overlap_file, "r") as f:
for line in f:
if line.strip() and not line.startswith("#"):
try:
index = int(line.strip().split()[0]) # Assume index is first column
indices_to_exclude.add(index)
except ValueError:
continue # Skip non-numeric or malformed lines
#Step 2: Read the stats file and extract header + table
with open(lh_stats_file, "r") as f:
lines = f.readlines()
#Find column header line
col_header_index = None
for i, line in enumerate(lines):
if line.startswith("# ColHeaders"):
col_header_index = i
break
if col_header_index is None:
raise ValueError("'# ColHeaders' not found in the stats file.")
# Get actual column names from header
columns = lines[col_header_index].strip().replace("# ColHeaders", "").split()
# Parse table lines into DataFrame
table_data = lines[col_header_index + 1:]
df = pd.read_csv(StringIO("".join(table_data)), delim_whitespace=True, names=columns)
# Step 3: Filter out rows with Index in exclusion list
df_filtered = df[~df["Index"].isin(indices_to_exclude)]
# Step 4: Write back the filtered stats file
with open(output_file, "w") as f:
# Write all lines before table
for line in lines[:col_header_index + 1]:
f.write(line)
# Write filtered table
df_filtered.to_csv(f, sep="\t", index=False, header=False)
print(f"Filtered stats saved to: {output_file}")
# %%
import pandas as pd
from io import StringIO
# Input files
lh_stats_file = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats/rh.w-g.pct.stats"
rh_overlap_file = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats/rh.lesion_roi_overlap_filtered.txt"
output_file = lh_stats_file.replace(".stats", ".filtered.stats")
# Step 1: Read the overlap file and collect Index values
indices_to_exclude = set()
with open(rh_overlap_file, "r") as f:
for line in f:
if line.strip() and not line.startswith("#"):
try:
index = int(line.strip().split()[0]) # Assume index is first column
indices_to_exclude.add(index)
except ValueError:
continue # Skip non-numeric or malformed lines
# Step 2: Read the stats file and extract header + table
with open(lh_stats_file, "r") as f:
lines = f.readlines()
# Find column header line
col_header_index = None
for i, line in enumerate(lines):
if line.startswith("# ColHeaders"):
col_header_index = i
break
if col_header_index is None:
raise ValueError("'# ColHeaders' not found in the stats file.")
# Get actual column names from header
columns = lines[col_header_index].strip().replace("# ColHeaders", "").split()
# Parse table lines into DataFrame
table_data = lines[col_header_index + 1:]
df = pd.read_csv(StringIO("".join(table_data)), delim_whitespace=True, names=columns)
# Step 3: Filter out rows with Index in exclusion list
df_filtered = df[~df["Index"].isin(indices_to_exclude)]
# Step 4: Write back the filtered stats file
with open(output_file, "w") as f:
# Write all lines before table
for line in lines[:col_header_index + 1]:
f.write(line)
# Write filtered table
df_filtered.to_csv(f, sep="\t", index=False, header=False)
print(f" Filtered stats saved to: {output_file}")
# %% [markdown]
# # Asymmetry Index
# %% [markdown]
# # FCD, Intensity
# %%
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# Define the main subject directory
main_folder = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats"
# Function to extract Mean intensity values from FreeSurfer stats files
def extract_intensity_values(stats_file):
extracted_values = {}
with open(stats_file, "r") as file:
for line in file:
if line.startswith("#") or not line.strip():
continue # Skip headers and blank lines
parts = line.split()
if len(parts) >= 6:
region = parts[4] # 5th column: StructName
try:
value = float(parts[5]) # 6th column: Mean intensity
extracted_values[region] = value
except ValueError:
extracted_values[region] = np.nan
return extracted_values
# Input paths for left and right hemisphere stats
lh_intensity_path = os.path.join(main_folder, "lh.w-g.pct.filtered.stats")
rh_intensity_path = os.path.join(main_folder, "rh.w-g.pct.filtered.stats")
# Extract data
lh_intensity_data = extract_intensity_values(lh_intensity_path)
rh_intensity_data = extract_intensity_values(rh_intensity_path)
# Compute AI for each ROI
roi_data = []
for region in set(lh_intensity_data) | set(rh_intensity_data):
left_value = lh_intensity_data.get(region, np.nan)
right_value = rh_intensity_data.get(region, np.nan)
if not np.isnan(left_value) and not np.isnan(right_value) and (left_value + right_value) != 0:
ai_value = (left_value - right_value) / (left_value + right_value)
else:
ai_value = np.nan
roi_data.append({
"StructName": region,
"Asymmetry_Index": ai_value,
"Left": left_value,
"Right": right_value
})
# Convert to DataFrame and save
df = pd.DataFrame(roi_data)
output_csv_path = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146_intensity_freesurfer_stats.csv"
df.to_csv(output_csv_path, index=False)
print(f"Intensity Asymmetry Index saved to: {output_csv_path}")
# Visualization
plt.figure(figsize=(20, 8))
sns.barplot(x="StructName", y="Asymmetry_Index", data=df.sort_values("Asymmetry_Index"), palette="coolwarm")
plt.title("Intensity Asymmetry Index (AI) per Brain Region in sub-00146", fontsize=18)
plt.xlabel("Brain Structure", fontsize=16)
plt.ylabel("Asymmetry Index (AI)", fontsize=16)
plt.xticks(rotation=90, fontsize=12)
plt.axhline(0, color='gray', linestyle='--')
plt.tight_layout()
plt.show()
# %% [markdown]
# # FCD, other structural measurements : (thickness, volume, curvature, surface area)
# %%
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# === CONFIGURATION ===
stats_dir = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna/sub-00146/stats_cleaned"
lh_stats = os.path.join(stats_dir, "lh.aparc.stats")
rh_stats = os.path.join(stats_dir, "rh.aparc.stats")
# Extract subject ID from path
subject_id = os.path.basename(os.path.dirname(stats_dir))
# Define metric columns (0-based)
metric_columns = {
"SurfArea": 2,
"GrayVol": 3,
"ThickAvg": 4,
"MeanCurv": 6
}
# === FUNCTIONS ===
def extract_metric(file_path, metric_name):
metric_idx = metric_columns[metric_name]
data = {}
with open(file_path, "r") as file:
for line in file:
if line.startswith("#") or not line.strip():
continue
parts = line.split()
roi = parts[0]
try:
value = float(parts[metric_idx])
data[roi] = value
except ValueError:
data[roi] = np.nan
return data
def plot_ai(metric_name):
lh_data = extract_metric(lh_stats, metric_name)
rh_data = extract_metric(rh_stats, metric_name)
roi_data = []
for roi in sorted(set(lh_data.keys()) | set(rh_data.keys())):
left = lh_data.get(roi, np.nan)
right = rh_data.get(roi, np.nan)
if not np.isnan(left) and not np.isnan(right) and (left + right) != 0:
ai = (left - right) / (left + right)
else:
ai = np.nan
roi_data.append({"StructName": roi, "Left": left, "Right": right, "Asymmetry_Index": ai})
df = pd.DataFrame(roi_data)
# Save to CSV with subject ID
output_file = f"{subject_id}_{metric_name}_asymmetry.csv"
output_path = os.path.join(stats_dir, output_file)
df.to_csv(output_path, index=False)
print(f"Saved {metric_name} AI to: {output_path}")
# Plot
plt.figure(figsize=(20, 8))
sns.barplot(x="StructName", y="Asymmetry_Index", data=df.sort_values("Asymmetry_Index"), palette="coolwarm")
plt.title(f"Asymmetry Index (AI) per Brain Region - {metric_name}", fontsize=18)
plt.xlabel("Brain Structure", fontsize=14)
plt.ylabel("Asymmetry Index (AI)", fontsize=14)
plt.xticks(rotation=90, fontsize=10)
plt.axhline(0, color='gray', linestyle='--')
plt.tight_layout()
plt.show()
# === Example run ===
for metric in metric_columns.keys():
plot_ai(metric)
# %% [markdown]
# # Healthy individuals: intensity
# %%
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# Define the main subject directory
main_folder = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_hc/sub-00170/stats"
# Function to extract Mean intensity values from FreeSurfer stats files
def extract_intensity_values(stats_file):
extracted_values = {}
with open(stats_file, "r") as file:
for line in file:
if line.startswith("#") or not line.strip():
continue # Skip headers and blank lines
parts = line.split()
if len(parts) >= 6:
region = parts[4] # 5th column: StructName
try:
value = float(parts[5]) # 6th column: Mean intensity
extracted_values[region] = value
except ValueError:
extracted_values[region] = np.nan
return extracted_values
# Input paths for left and right hemisphere stats
lh_intensity_path = os.path.join(main_folder, "lh.w-g.pct.stats")
rh_intensity_path = os.path.join(main_folder, "rh.w-g.pct.stats")
# Extract data
lh_intensity_data = extract_intensity_values(lh_intensity_path)
rh_intensity_data = extract_intensity_values(rh_intensity_path)
# Compute asymmetry index for each ROI
roi_data = []
for region in set(lh_intensity_data) | set(rh_intensity_data):
left_value = lh_intensity_data.get(region, np.nan)
right_value = rh_intensity_data.get(region, np.nan)
if not np.isnan(left_value) and not np.isnan(right_value) and (left_value + right_value) != 0:
ai_value =(left_value - right_value) / (left_value + right_value)
else:
ai_value = np.nan
roi_data.append({
"StructName": region,
"Asymmetry_Index": ai_value,
"Left": left_value,
"Right": right_value
})
# Convert to DataFrame and save
df = pd.DataFrame(roi_data)
output_csv_path = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_hc/sub-00170_intensity_freesurfer_stats.csv"
df.to_csv(output_csv_path, index=False)
print(f"Intensity Asymmetry Index saved to: {output_csv_path}")
# Visualization
plt.figure(figsize=(20, 8))
sns.barplot(x="StructName", y="Asymmetry_Index", data=df.sort_values("Asymmetry_Index"), palette="coolwarm")
plt.title("Intensity Asymmetry Index (AI) per Brain Region in sub-00170", fontsize=18)
plt.xlabel("Brain Structure", fontsize=16)
plt.ylabel("Asymmetry Index (AI)", fontsize=16)
plt.xticks(rotation=90, fontsize=12)
plt.axhline(0, color='gray', linestyle='--')
plt.tight_layout()
plt.show()
# %% [markdown]
# # Healthy individuals: other cortical measurements: (thickness, volume, curvature, surface area)
# %%
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# === CONFIGURATION ===
stats_dir = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_hc/sub-00170/stats"
lh_stats = os.path.join(stats_dir, "lh.aparc.stats")
rh_stats = os.path.join(stats_dir, "rh.aparc.stats")
# Extract subject ID from path
subject_id = os.path.basename(os.path.dirname(stats_dir))
# Define metric columns (0-based)
metric_columns = {
"SurfArea": 2,
"GrayVol": 3,
"ThickAvg": 4,
"MeanCurv": 6
}
# === FUNCTIONS ===
def extract_metric(file_path, metric_name):
metric_idx = metric_columns[metric_name]
data = {}
with open(file_path, "r") as file:
for line in file:
if line.startswith("#") or not line.strip():
continue
parts = line.split()
roi = parts[0]
try:
value = float(parts[metric_idx])
data[roi] = value
except ValueError:
data[roi] = np.nan
return data
def plot_ai(metric_name):
lh_data = extract_metric(lh_stats, metric_name)
rh_data = extract_metric(rh_stats, metric_name)
roi_data = []
for roi in sorted(set(lh_data.keys()) | set(rh_data.keys())):
left = lh_data.get(roi, np.nan)
right = rh_data.get(roi, np.nan)
if not np.isnan(left) and not np.isnan(right) and (left + right) != 0:
ai = (left - right) / (left + right)
else:
ai = np.nan
roi_data.append({"StructName": roi, "Left": left, "Right": right, "Asymmetry_Index": ai})
df = pd.DataFrame(roi_data)
# Save to CSV with subject ID
output_file = f"{subject_id}_{metric_name}_asymmetry.csv"
output_path = os.path.join(stats_dir, output_file)
df.to_csv(output_path, index=False)
print(f"Saved {metric_name} AI to: {output_path}")
# Plot
plt.figure(figsize=(20, 8))
sns.barplot(x="StructName", y="Asymmetry_Index", data=df.sort_values("Asymmetry_Index"), palette="coolwarm")
plt.title(f"Asymmetry Index (AI) per Brain Region - {metric_name}", fontsize=18)
plt.xlabel("Brain Structure", fontsize=14)
plt.ylabel("Asymmetry Index (AI)", fontsize=14)
plt.xticks(rotation=90, fontsize=10)
plt.axhline(0, color='gray', linestyle='--')
plt.tight_layout()
plt.show()
# === Example run ===
for metric in metric_columns.keys():
plot_ai(metric)
# %% [markdown]
# # Independence check
# %%
import os
import pandas as pd
import numpy as np
from scipy.stats import spearmanr, chi2_contingency
import warnings
warnings.filterwarnings("ignore")
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/participants.aparc.tsv")
# === Load participant metadata ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# --- Clean metadata ---
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['group'] = participants_df['group'].astype(str).str.lower()
participants_df['histopathology'] = participants_df['histopathology'].astype(str).str.lower()
participants_df.dropna(subset=['sex', 'group', 'age_scan'], inplace=True)
# Encode sex for numeric correlations
sex_map = {'F': 0, 'M': 1}
participants_df['sex_numeric'] = participants_df['sex'].map(sex_map)
# === Prepare results container ===
results = []
# === Robust independence check function ===
def check_independence(df, numeric_col=None, cat_col=None, contrast_name=""):
n_total = len(df)
n_per_group = df[cat_col].value_counts().to_dict() if cat_col else None
if df.empty:
results.append({
"contrast": contrast_name,
"test_type": "info",
"statistic": None,
"p_value": None,
"note": "No data available"
})
return
# Spearman correlation
if numeric_col is not None and 'age_scan' in df.columns:
if df[numeric_col].nunique() > 1:
corr, pval = spearmanr(df['age_scan'], df[numeric_col])
results.append({
"contrast": contrast_name,
"test_type": "spearman",
"statistic": corr,
"p_value": pval,
"note": f"n_total={n_total}, n_per_group={n_per_group}"
})
else:
results.append({
"contrast": contrast_name,
"test_type": "spearman",
"statistic": None,
"p_value": None,
"note": f"Only one unique value, n_total={n_total}, n_per_group={n_per_group}"
})
# Chi-square test
if cat_col is not None:
contingency = pd.crosstab(df['sex'], df[cat_col])
if contingency.size == 0 or contingency.shape[0] < 2 or contingency.shape[1] < 2:
results.append({
"contrast": contrast_name,
"test_type": "chi2",
"statistic": None,
"p_value": None,
"note": f"Not enough data, n_total={n_total}, n_per_group={n_per_group}"
})
else:
chi2, p, _, _ = chi2_contingency(contingency)
results.append({
"contrast": contrast_name,
"test_type": "chi2",
"statistic": chi2,
"p_value": p,
"note": f"n_total={n_total}, n_per_group={n_per_group}"
})
# === 1. HC vs FCD (use 'group' column) ===
hc_vs_fcd = participants_df[participants_df['group'].isin(['hc','fcd'])].copy()
hc_vs_fcd['group_numeric'] = hc_vs_fcd['group'].map({'hc':0,'fcd':1})
check_independence(hc_vs_fcd, numeric_col='group_numeric', cat_col='group', contrast_name="HC vs FCD")
# === 2-4. FCD subtype contrasts (use 'histopathology' column) ===
fcd_subtypes = participants_df[participants_df['histopathology'].isin(['fcdna','iia','iib'])].copy()
subtype_pairs = [('fcdna','iia'), ('fcdna','iib'), ('iia','iib')]
for sub1, sub2 in subtype_pairs:
df_pair = fcd_subtypes[fcd_subtypes['histopathology'].isin([sub1, sub2])].copy()
df_pair['histopathology_numeric'] = df_pair['histopathology'].map({sub1:0, sub2:1})
check_independence(df_pair, numeric_col='histopathology_numeric', cat_col='histopathology',
contrast_name=f"{sub1.upper()} vs {sub2.upper()}")
# === Save results to CSV and Excel ===
results_df = pd.DataFrame(results)
csv_file = os.path.join(data_root, "independence_results.csv")
excel_file = os.path.join(data_root, "independence_results.xlsx")
results_df.to_csv(csv_file, index=False)
results_df.to_excel(excel_file, index=False)
print(f"Independence check results saved to:\nCSV: {csv_file}\nExcel: {excel_file}")
# === Show results after running ===
print("\n=== Independence Check Results ===")
print(results_df)
# %% [markdown]
# ## Beta regression for intensity FCD vs HC
# %%
pip install rpy2
# %%
import os
import pandas as pd
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
import warnings
warnings.filterwarnings("ignore")
# ==== Load R package betareg ====
utils = importr('utils')
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# Clean metadata
participants_df['group'] = participants_df['group'].replace({'fcdlla': 'fcd', 'fcdllb': 'fcd', 'fcdna': 'fcd'})
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['group'] = participants_df['group'].astype(str).str.lower()
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
group_map = {'hc': 0, 'fcd': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
participants_df['group'] = participants_df['group'].map(group_map)
participants_df.dropna(subset=['sex', 'group', 'age_scan'], inplace=True)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and group:", spearmanr(participants_df['age_scan'], participants_df['group']))
contingency_sex_group = pd.crosstab(participants_df['sex'], participants_df['group'])
chi2, p, _, _ = chi2_contingency(contingency_sex_group)
print("Chi-square test between sex and group: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna", "asymmetry_hc"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
continue
for file_name in os.listdir(group_path):
if not file_name.endswith("_intensity_freesurfer_stats.csv"):
continue
file_path = os.path.join(group_path, file_name)
subject = file_name.replace("_intensity_freesurfer_stats.csv", "")
try:
df = pd.read_csv(file_path)
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'group', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Intensity Asymmetry FCD vs HC")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %%
!pip install nilearn
# %% [markdown]
# ## Beta regression for thickness FCD vs HC
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# Clean metadata
participants_df['group'] = participants_df['group'].replace({'fcdlla': 'fcd', 'fcdllb': 'fcd', 'fcdna': 'fcd'})
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['group'] = participants_df['group'].astype(str).str.lower()
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
group_map = {'hc': 0, 'fcd': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
participants_df['group'] = participants_df['group'].map(group_map)
participants_df.dropna(subset=['sex', 'group', 'age_scan'], inplace=True)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and group:", spearmanr(participants_df['age_scan'], participants_df['group']))
contingency_sex_group = pd.crosstab(participants_df['sex'], participants_df['group'])
chi2, p, _, _ = chi2_contingency(contingency_sex_group)
print("Chi-square test between sex and group: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna", "asymmetry_hc"]
all_data = []
# Load thickness asymmetry data from nested structure
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
if group_folder == "asymmetry_hc":
stats_folder = os.path.join(subj_folder_path, "stats")
else:
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_ThickAvg_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'group', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Thickness Asymmetry FCD vs HC")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for volume FCD vs HC
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# Clean metadata
participants_df['group'] = participants_df['group'].replace({'fcdlla': 'fcd', 'fcdllb': 'fcd', 'fcdna': 'fcd'})
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['group'] = participants_df['group'].astype(str).str.lower()
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
group_map = {'hc': 0, 'fcd': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
participants_df['group'] = participants_df['group'].map(group_map)
participants_df.dropna(subset=['sex', 'group', 'age_scan'], inplace=True)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and group:", spearmanr(participants_df['age_scan'], participants_df['group']))
contingency_sex_group = pd.crosstab(participants_df['sex'], participants_df['group'])
chi2, p, _, _ = chi2_contingency(contingency_sex_group)
print("Chi-square test between sex and group: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna", "asymmetry_hc"]
all_data = []
# Load volume asymmetry data from nested structure
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
if group_folder == "asymmetry_hc":
stats_folder = os.path.join(subj_folder_path, "stats")
else:
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_GrayVol_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'group', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Volume Asymmetry FCD vs HC")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for curvature FCD vs HC
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# Clean metadata
participants_df['group'] = participants_df['group'].replace({'fcdlla': 'fcd', 'fcdllb': 'fcd', 'fcdna': 'fcd'})
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['group'] = participants_df['group'].astype(str).str.lower()
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
group_map = {'hc': 0, 'fcd': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
participants_df['group'] = participants_df['group'].map(group_map)
participants_df.dropna(subset=['sex', 'group', 'age_scan'], inplace=True)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and group:", spearmanr(participants_df['age_scan'], participants_df['group']))
contingency_sex_group = pd.crosstab(participants_df['sex'], participants_df['group'])
chi2, p, _, _ = chi2_contingency(contingency_sex_group)
print("Chi-square test between sex and group: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna", "asymmetry_hc"]
all_data = []
# Load curvature asymmetry data from nested structure
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
if group_folder == "asymmetry_hc":
stats_folder = os.path.join(subj_folder_path, "stats")
else:
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_MeanCurv_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'group', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Curvature Asymmetry FCD vs HC")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for surface area FCD vs HC
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# Clean metadata
participants_df['group'] = participants_df['group'].replace({'fcdlla': 'fcd', 'fcdllb': 'fcd', 'fcdna': 'fcd'})
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['group'] = participants_df['group'].astype(str).str.lower()
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
group_map = {'hc': 0, 'fcd': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
participants_df['group'] = participants_df['group'].map(group_map)
participants_df.dropna(subset=['sex', 'group', 'age_scan'], inplace=True)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and group:", spearmanr(participants_df['age_scan'], participants_df['group']))
contingency_sex_group = pd.crosstab(participants_df['sex'], participants_df['group'])
chi2, p, _, _ = chi2_contingency(contingency_sex_group)
print("Chi-square test between sex and group: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna", "asymmetry_hc"]
all_data = []
# Load surface area asymmetry data from nested structure
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
if group_folder == "asymmetry_hc":
stats_folder = os.path.join(subj_folder_path, "stats")
else:
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_SurfArea_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'group', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Surface Area Asymmetry FCD vs HC")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for intensity FCD subtypes
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
continue
for file_name in os.listdir(group_path):
if not file_name.endswith("_intensity_freesurfer_stats.csv"):
continue
file_path = os.path.join(group_path, file_name)
subject = file_name.replace("_intensity_freesurfer_stats.csv", "")
try:
df = pd.read_csv(file_path)
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Intensity Asymmetry between FCD subtypes (FCDIIa, FCDIIb, and non-operated FCD)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for thickness FCD subtypes
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_ThickAvg_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Thickness Asymmetry between FCD subtypes (FCDIIa, FCDIIb, and non-operated FCD)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for volume FCD subtypes
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_GrayVol_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Volume Asymmetry between FCD subtypes (FCDIIa, FCDIIb, and non-operated FCD)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for curvature FCD subtypes
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_MeanCurv_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Curvature Asymmetry between FCD subtypes (FCDIIa, FCDIIb, and non-operated FCD)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for surface area FCD subtypes
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_SurfArea_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Surface Area Asymmetry between FCD subtypes (FCDIIa, FCDIIb, and non-operated FCD)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for intensity FCDIIb vs FCDIIa
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
continue
for file_name in os.listdir(group_path):
if not file_name.endswith("_intensity_freesurfer_stats.csv"):
continue
file_path = os.path.join(group_path, file_name)
subject = file_name.replace("_intensity_freesurfer_stats.csv", "")
try:
df = pd.read_csv(file_path)
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Intensity Asymmetry between FCD subtypes (FCDIIa, FCDIIb)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for thickness FCDIIb vs FCDIIa
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_ThickAvg_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Thickness Asymmetry between FCD subtypes (FCDIIa, FCDIIb)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for volume FCDIIb vs FCDIIa
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_GrayVol_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Volume Asymmetry between FCD subtypes (FCDIIa, FCDIIb)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for curvature FCDIIb vs FCDIIa
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_MeanCurv_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Curvature Asymmetry between FCD subtypes (FCDIIa, FCDIIb)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %% [markdown]
# ## Beta regression for surface area FCDIIb vs FCDIIa
# %%
import os
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import spearmanr, chi2_contingency, kstest
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multitest import multipletests
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import r, globalenv, Formula
from rpy2.robjects.packages import importr
betareg = importr('betareg')
import warnings
warnings.filterwarnings("ignore")
# ==== 3. Load R packages ====
utils = importr('utils')
utils.install_packages('betareg') # will skip if already installed
betareg = importr('betareg')
# === Base directory and metadata path ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
metadata_csv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv")
metadata_tsv = os.path.join(data_root, "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv")
# === Load participant metadata safely ===
try:
if os.path.exists(metadata_csv):
participants_df = pd.read_csv(metadata_csv, sep=None, engine='python', on_bad_lines='skip')
elif os.path.exists(metadata_tsv):
participants_df = pd.read_csv(metadata_tsv, sep='\t', engine='python', on_bad_lines='skip')
else:
raise FileNotFoundError("Neither CSV nor TSV metadata file found.")
except Exception as e:
raise RuntimeError(f"Failed to load metadata: {e}")
# === Standardize and encode columns ===
# Drop rows missing important fields first
participants_df.dropna(subset=['sex', 'age_scan','histopathology'], inplace=True)
# Clean metadata
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# Encode categorical variables
sex_map = {'F': 0, 'M': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
# Encode histopathology (use instead of group)
participants_df['histopathology'] = participants_df['histopathology'].str.lower().map({'iia': 0, 'iib': 1, 'na': 2})
# Ensure 'participant_id' exists and is string
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
# === Independence Checks ===
print("\n=== Independence Checks ===")
print("Spearman correlation between age_scan and sex:", spearmanr(participants_df['age_scan'], participants_df['sex']))
print("Spearman correlation between age_scan and histopathology:", spearmanr(participants_df['age_scan'], participants_df['histopathology']))
contingency_sex_histopathology = pd.crosstab(participants_df['sex'], participants_df['histopathology'])
chi2, p, _, _ = chi2_contingency(contingency_sex_histopathology)
print("Chi-square test between sex and histopathology: chi2 =", chi2, ", p =", p)
# === Load asymmetry data ===
group_folders = ["asymmetry_fcdlla", "asymmetry_fcdllb"]
all_data = []
for group_folder in group_folders:
group_path = os.path.join(data_root, group_folder)
if not os.path.exists(group_path):
print(f"Missing: {group_path}")
continue
# Iterate over subject folders inside group folder
for subject in os.listdir(group_path):
subj_folder_path = os.path.join(group_path, subject)
if not os.path.isdir(subj_folder_path):
continue
# Use 'stats_cleaned' except for HC group which uses 'stats'
stats_folder = os.path.join(subj_folder_path, "stats_cleaned")
file_name = f"{subject}_SurfArea_asymmetry.csv"
file_path = os.path.join(stats_folder, file_name)
if not os.path.exists(file_path):
print(f"Not found: {file_path}")
continue
try:
df = pd.read_csv(file_path)
# Check necessary columns
if not all(col in df.columns for col in ['StructName', 'Asymmetry_Index']):
print(f"Missing columns in {file_path}")
continue
df['participant_id'] = subject
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
all_data.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
except Exception as e:
print(f"Failed loading {file_path}: {e}")
# Merge all data
asymmetry_df = pd.concat(all_data, ignore_index=True)
full_df = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
full_df.dropna(subset=['Unsigned_Asymmetry', 'sex', 'histopathology', 'age_scan'], inplace=True)
# === Check response distribution and transform if needed ===
plt.figure(figsize=(8, 4))
sns.histplot(full_df['Unsigned_Asymmetry'], kde=True)
plt.title("Distribution of Surface Area Asymmetry between FCD subtypes (FCDIIa, FCDIIb)")
plt.show()
ks_stat, ks_p = kstest(full_df['Unsigned_Asymmetry'], 'norm')
if ks_p < 0.05:
full_df['Unsigned_Asymmetry'] = np.sqrt(full_df['Unsigned_Asymmetry'])
print("Data was not normal. Applied square root transformation.")
# %%
import rpy2.robjects.packages as rpackages
utils = rpackages.importr('utils')
utils.install_packages('statmod')
# %%
import rpy2.robjects as robjects
robjects.r('install.packages("numDeriv", repos="https://cloud.r-project.org/")')
robjects.r('install.packages("betareg", repos="https://cloud.r-project.org/")')
# %% [markdown]
# # FDR for all predictors * measures* ROIs, FCD vs HC Storey’s q-values FDR Correction
# %%
from rpy2.robjects.packages import importr
utils = importr("utils")
utils.install_packages("qvalue")
# %% [markdown]
# # For group only FCD vs HC
# %%
import os
from pathlib import Path
import logging
import pandas as pd
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from statsmodels.stats.multitest import multipletests
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore")
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri, r, Formula, globalenv
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
import rpy2.robjects.packages as rpackages
from rpy2.robjects.packages import importr
utils = importr('utils')
utils.chooseCRANmirror(ind=1)
# Try BiocManager if qvalue is not on CRAN
if not rpackages.isinstalled('qvalue'):
try:
if not rpackages.isinstalled('BiocManager'):
utils.install_packages('BiocManager')
biocmanager = importr('BiocManager')
biocmanager.install('qvalue')
except Exception as e:
raise RuntimeError(f"Failed to install qvalue via BiocManager: {e}")
# Now import
qvalue_pkg = importr('qvalue')
utils = importr('utils')
utils.chooseCRANmirror(ind=1)
if not rpackages.isinstalled('qvalue'):
try:
utils.install_packages('qvalue', repos='https://cloud.r-project.org', Ncpus=4)
except Exception as e:
raise RuntimeError(f"Failed to install qvalue via CRAN: {e}")
# Logging
logging.basicConfig(level=logging.INFO, format="%(asctime)s [%(levelname)s] %(message)s")
logger = logging.getLogger("betareg_global")
# R packages
utils = importr('utils')
utils.install_packages('betareg')
try:
betareg = importr('betareg')
except Exception as e:
raise RuntimeError("R package 'betareg' not available. Please install it in R.") from e
# User paths
DATA_ROOT = Path("/Volumes/groups/tohkagroup/Bonn_Epilepsy")
METADATA_XLSX = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.xlsx"
METADATA_CSV = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv"
METADATA_TSV = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv"
OUT_ROOT = DATA_ROOT / "betareg_global_results_fcd_vs_hc_storey_qvalue_&_benjamini_hochberg_ group_only"
OUT_ROOT.mkdir(parents=True, exist_ok=True)
# Load participant metadata
participants_df = None
if METADATA_XLSX.exists():
try:
participants_df = pd.read_excel(METADATA_XLSX, engine='openpyxl')
logger.info(f"Loaded metadata from Excel: {METADATA_XLSX}")
except Exception as e:
logger.warning(f"Failed Excel read {METADATA_XLSX}: {e}")
if participants_df is None and METADATA_CSV.exists():
try:
participants_df = pd.read_csv(METADATA_CSV, sep=';', engine='python', on_bad_lines='skip')
logger.info(f"Loaded metadata from CSV: {METADATA_CSV}")
except Exception as e:
logger.warning(f"Failed CSV read {METADATA_CSV}: {e}")
if participants_df is None and METADATA_TSV.exists():
try:
participants_df = pd.read_csv(METADATA_TSV, sep='\t', engine='python', on_bad_lines='skip')
logger.info(f"Loaded metadata from TSV: {METADATA_TSV}")
except Exception as e:
logger.warning(f"Failed TSV read {METADATA_TSV}: {e}")
if participants_df is None:
raise FileNotFoundError("Metadata not found or unreadable.")
# Clean metadata
if 'participant_id' not in participants_df.columns and 'subject_id' in participants_df.columns:
participants_df = participants_df.rename(columns={'subject_id': 'participant_id'})
participants_df['group'] = participants_df['group'].astype(str).str.lower().replace({
'fcdlla': 'fcd', 'fcdllb': 'fcd', 'fcdna': 'fcd'
})
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
sex_map = {'F': 0, 'M': 1}
group_map = {'hc': 0, 'fcd': 1}
participants_df['sex'] = participants_df['sex'].map(sex_map)
participants_df['group'] = participants_df['group'].map(group_map)
participants_df = participants_df.dropna(subset=['participant_id', 'sex', 'group', 'age_scan'])
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
logger.info(f"Loaded metadata: {len(participants_df)} participants after cleaning.")
# Settings
MEASURES = ["intensity", "thickness", "volume", "curvature", "surface_area"]
MEASURE_FILE_PATTERNS = {
"intensity": "_intensity_freesurfer_stats.csv",
"thickness": "_ThickAvg_asymmetry.csv",
"volume": "_GrayVol_asymmetry.csv",
"curvature": "_MeanCurv_asymmetry.csv",
"surface_area": "_SurfArea_asymmetry.csv",
}
GROUP_FOLDERS = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna", "asymmetry_hc"]
MIN_SUBJ_PER_ROI = 12
EPS = 1e-4
ALPHA = 0.05
CI_Z = 1.96
# ---------- Helper functions ----------
def collect_intensity_files():
files = []
for folder in GROUP_FOLDERS:
gp = DATA_ROOT / folder
if not gp.exists():
continue
for fn in os.listdir(gp):
if fn.endswith(MEASURE_FILE_PATTERNS['intensity']):
files.append(gp / fn)
return files
def collect_nested_measure_files(measure_key):
files = []
pattern_suffix = MEASURE_FILE_PATTERNS[measure_key]
for folder in GROUP_FOLDERS:
gp = DATA_ROOT / folder
if not gp.exists():
continue
for subj in os.listdir(gp):
subj_folder = gp / subj
if not subj_folder.is_dir():
continue
stats_folder = subj_folder / ("stats" if folder=="asymmetry_hc" else "stats_cleaned")
if not stats_folder.exists():
stats_folder = subj_folder / "stats"
if not stats_folder.exists():
continue
fpath = stats_folder / f"{subj}{pattern_suffix}"
if fpath.exists():
files.append(fpath)
return files
# ---------- Main loop ----------
all_results = []
failures = []
for measure in MEASURES:
logger.info(f"Processing measure: {measure}")
measure_out = OUT_ROOT / measure
measure_out.mkdir(parents=True, exist_ok=True)
files = collect_intensity_files() if measure=="intensity" else collect_nested_measure_files(measure)
if not files:
logger.warning(f"No files found for {measure}")
continue
rows = []
for fpath in files:
try:
df = pd.read_csv(fpath)
except:
continue
asym_col = 'Asymmetry_Index' if 'Asymmetry_Index' in df.columns else ('AsymmetryIndex' if 'AsymmetryIndex' in df.columns else None)
if asym_col is None:
continue
fname = fpath.name
participant_id = fname.replace(MEASURE_FILE_PATTERNS[measure], '') if measure=='intensity' else fname[:-len(MEASURE_FILE_PATTERNS[measure])]
df = df[['StructName', asym_col]].rename(columns={asym_col:'Asymmetry_Index'})
df['participant_id'] = participant_id
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
rows.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
if not rows:
continue
asymmetry_df = pd.concat(rows, ignore_index=True)
merged = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
merged = merged.dropna(subset=['Unsigned_Asymmetry','sex','group','age_scan'])
if merged.empty:
continue
# Save distribution plot
plt.figure(figsize=(6,3))
sns.histplot(merged['Unsigned_Asymmetry'], bins=40)
plt.title(f"{measure} - Unsigned_Asymmetry distribution")
plt.tight_layout()
plt.savefig(measure_out / f"{measure}_scaled_dist.png")
plt.close()
roi_list = merged['StructName'].unique()
formula = Formula('Unsigned_Asymmetry ~ group + sex + age_scan')
for roi in tqdm(roi_list, desc=f"{measure} ROIs"):
roi_df = merged[merged['StructName']==roi].copy()
nsubj = roi_df['participant_id'].nunique()
if nsubj < MIN_SUBJ_PER_ROI or roi_df['Unsigned_Asymmetry'].nunique()<3:
failures.append({'Measure':measure,'ROI':roi,'Reason':'too_few_or_low_variability','N':nsubj})
continue
try:
with localconverter(ro.default_converter + pandas2ri.converter):
rdf = ro.conversion.py2rpy(roi_df[['Unsigned_Asymmetry','group','sex','age_scan']])
globalenv['rdf'] = rdf
model = betareg.betareg(formula, data=rdf)
summary_model = r.summary(model)
# --- Residual plots ---
try:
r_resid = r.residuals(model)
with localconverter(ro.default_converter + pandas2ri.converter):
py_resid = ro.conversion.rpy2py(r_resid)
resid_series = pd.Series(py_resid).astype(float)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Histogram
axes[0].hist(resid_series, bins=30, color='skyblue', edgecolor='black')
axes[0].set_title(f"{measure} {roi} residuals histogram")
# QQ plot
import statsmodels.api as sm
sm.graphics.qqplot(resid_series, line='45', ax=axes[1])
axes[1].set_title(f"{measure} {roi} residuals QQ")
plt.tight_layout()
plt.savefig(measure_out / f"{roi}_residuals.png", dpi=200)
plt.close()
except Exception as e:
logger.warning(f"Residuals plot failed for {measure} {roi}: {e}")
coef_table = summary_model.rx2('coefficients').rx2('mean')
coef_names = list(coef_table.rownames)
for name in coef_names:
est = float(coef_table.rx(name, 'Estimate')[0])
se = float(coef_table.rx(name, 'Std. Error')[0])
z = float(coef_table.rx(name, 'z value')[0])
p = float(coef_table.rx(name, 'Pr(>|z|)')[0])
ci_low = est - CI_Z * se
ci_high = est + CI_Z * se
all_results.append({
'Measure': measure,
'ROI': roi,
'Predictor': name,
'Estimate': est,
'StdErr': se,
'z': z,
'Pval': p, # raw p-value
'CI_low': ci_low,
'CI_high': ci_high,
'Pseudo_R2': float(summary_model.rx2('pseudo.r.squared')[0]) if 'pseudo.r.squared' in summary_model.names else np.nan,
'N_subjects': nsubj
})
except Exception as e:
failures.append({'Measure':measure,'ROI':roi,'Reason':f"betareg_failed:{e}",'N':nsubj})
continue
# ---------- Results and corrections ----------
results_df = pd.DataFrame(all_results)
fails_df = pd.DataFrame(failures)
if results_df.empty:
raise RuntimeError("No model results collected.")
# --- Keep raw p-values separately for both stages ---
results_df['Pval_raw_global'] = results_df['Pval'].astype(float)
results_df['Pval_raw_predictor'] = results_df['Pval'].astype(float)
# --- Global FDR (all ROIs and measures) but only for group/histopathology ---
mask_interest = results_df['Predictor'].isin(['group', 'histopathology'])
pvals_all = results_df.loc[mask_interest, 'Pval_raw_global'].astype(float).values
# Benjamini-Hochberg (FDR)
fdr_adj = multipletests(pvals_all, alpha=ALPHA, method='fdr_bh')[1]
results_df.loc[mask_interest, 'FDR_p_global'] = fdr_adj
# Bonferroni correction
bonf_adj = multipletests(pvals_all, alpha=ALPHA, method='bonferroni')[1]
results_df.loc[mask_interest, 'Bonferroni_p'] = bonf_adj
# Flags
results_df.loc[mask_interest, 'significant_FDR_global'] = fdr_adj < ALPHA
results_df.loc[mask_interest, 'significant_Bonferroni'] = bonf_adj < ALPHA
# --- Storey q-values (using R's qvalue package) for global correction ---
try:
qvalue_pkg = importr('qvalue')
with localconverter(ro.default_converter + pandas2ri.converter):
r_pvals = ro.FloatVector(pvals_all)
qobj = qvalue_pkg.qvalue(r_pvals)
qvalues = np.array(qobj.rx2('qvalues'))
results_df.loc[mask_interest, 'qvalue_global'] = qvalues
results_df.loc[mask_interest, 'significant_qvalue'] = qvalues < ALPHA
logger.info("Storey q-values successfully computed for group/histopathology only (global correction).")
except Exception as e:
logger.warning(f"Storey q-value computation failed for global correction: {e}")
results_df['qvalue_global'] = np.nan
results_df['significant_qvalue'] = False
logger.info("Storey q-values successfully computed and added.")
except Exception as e:
logger.warning(f"Storey q-value computation failed: {e}")
results_df['qvalue_global'] = np.nan
results_df['significant_qvalue'] = False
# --- Per-predictor FDR (only for 'group' or 'histopathology') ---
results_df['FDR_p_within_predictor'] = np.nan
results_df['qvalue_within_predictor'] = np.nan
results_df['significant_qvalue_within_predictor'] = False
predictors_to_correct = ['group', 'histopathology'] # <only group or histopathology
for (measure, predictor), grp in results_df.groupby(['Measure','Predictor']):
if predictor not in predictors_to_correct:
continue # skip all other predictors like age, sex, etc.
# FDR correction (Benjamini-Hochberg)
adj_bh = multipletests(grp['Pval_raw_predictor'].values, method='fdr_bh')[1]
results_df.loc[grp.index, 'FDR_p_within_predictor'] = adj_bh
# Storey q-value
try:
with localconverter(ro.default_converter + pandas2ri.converter):
r_pvals = ro.FloatVector(grp['Pval_raw_predictor'].values)
qobj = qvalue_pkg.qvalue(r_pvals)
qvals = np.array(qobj.rx2('qvalues'))
results_df.loc[grp.index, 'qvalue_within_predictor'] = qvals
results_df.loc[grp.index, 'significant_qvalue_within_predictor'] = qvals < ALPHA
except Exception as e:
logger.warning(f"Storey q-value failed for predictor {predictor} in {measure}: {e}")
# --- Save outputs ---
results_df.to_csv(OUT_ROOT / "all_measures_betareg_results_with_global_corrections.csv", index=False)
try:
results_df.to_excel(OUT_ROOT / "all_measures_betareg_results_with_global_corrections.xlsx", index=False)
except:
logger.warning("Could not write Excel (check openpyxl)")
fails_df.to_csv(OUT_ROOT / "model_failures.csv", index=False)
logger.info(f"Saved results and failures in {OUT_ROOT}")
# %% [markdown]
# # Plot for FCD vs HC only group
# %%
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from nilearn import datasets, plotting
from nibabel.freesurfer.io import read_annot
# === Directories ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/betareg_global_results_fcd_vs_hc_storey_qvalue_&_benjamini_hochberg_ group_only"
results_csv = os.path.join(data_root, "all_measures_betareg_results_with_global_corrections.csv")
output_dir = os.path.join(data_root, "brain_maps_hc_fcd_zscore_global_by_measure")
os.makedirs(output_dir, exist_ok=True)
# === Load results ===
df = pd.read_csv(results_csv)
# === Filter significant ROIs (global q < 0.05, exclude intercepts) ===
sig_df = df[
(df['qvalue_global'] < 0.05) &
(~df['Predictor'].str.lower().str.contains("intercept"))
]
if sig_df.empty:
print(" No significant ROIs found with qvalue_global < 0.05 — plotting all instead.")
sig_df = df.copy()
# === Load fsaverage5 + aparc labels ===
fsaverage = datasets.fetch_surf_fsaverage('fsaverage5')
labels_left, ctab_left, names_left = read_annot(os.path.join(data_root, "lh.aparc.annot"))
labels_right, ctab_right, names_right = read_annot(os.path.join(data_root, "rh.aparc.annot"))
names_left = [n.decode("utf-8").lower() for n in names_left]
names_right = [n.decode("utf-8").lower() for n in names_right]
# === Map ROI z-scores to surface ===
def roi_to_surface(df_pred):
roi_dict = {roi.lower(): z for roi, z in zip(df_pred['ROI'], df_pred['z'])}
data_left = np.zeros(len(labels_left))
data_right = np.zeros(len(labels_right))
for roi_name, z in roi_dict.items():
if roi_name in names_left:
idx = names_left.index(roi_name)
data_left[labels_left == idx] = z
if roi_name in names_right:
idx = names_right.index(roi_name)
data_right[labels_right == idx] = z
return data_left, data_right
# === Plotting function ===
def save_views(data, hemi, predictor, measure):
for view in ['lateral', 'medial']:
fname = os.path.join(output_dir, f"{measure}_{predictor}_{hemi}_{view}.png")
fig = plt.figure(facecolor='black')
ax = fig.add_subplot(111, projection='3d')
plotting.plot_surf_stat_map(
fsaverage[f'infl_{hemi}'], data,
hemi=hemi,
bg_map=fsaverage[f'sulc_{hemi}'],
cmap='bwr', colorbar=True,
bg_on_data=True, darkness=0.5, alpha=1.0, axes=ax,
vmin=-5, vmax=5, view=view
)
for child_ax in fig.get_axes():
child_ax.set_facecolor("black")
fig.savefig(fname, dpi=300, facecolor='black', bbox_inches='tight')
plt.close(fig)
# === Generate plots per measure ===
measures = ['intensity', 'thickness', 'volume', 'curvature', 'surface_area']
for measure in measures:
df_measure = sig_df[sig_df['Measure'].str.lower() == measure.lower()]
if df_measure.empty:
print(f" No significant results found for {measure}, skipping.")
continue
predictors = df_measure['Predictor'].unique()
# Step 1: generate surface images
for predictor in predictors:
df_pred = df_measure[df_measure['Predictor'] == predictor]
data_left, data_right = roi_to_surface(df_pred)
save_views(data_left, 'left', predictor, measure)
save_views(data_right, 'right', predictor, measure)
print("Saved all individual predictor brain maps to:", output_dir)
# Step 2: assemble composite for this measure
n_preds = len(predictors)
fig, axes = plt.subplots(n_preds, 4, figsize=(16, 4 * n_preds), facecolor='black')
if n_preds == 1:
axes = np.array([axes])
for i, predictor in enumerate(predictors):
for j, hemi in enumerate(['left', 'right']):
for k, view in enumerate(['lateral', 'medial']):
img_path = os.path.join(output_dir, f"{measure}_{predictor}_{hemi}_{view}.png")
if os.path.exists(img_path):
img = plt.imread(img_path)
axes[i, j*2 + k].imshow(img)
axes[i, j*2 + k].axis("off")
axes[i, j*2 + k].set_facecolor('black')
axes[i, j*2 + k].set_title(f"{predictor} ({hemi} {view})", color='yellow', fontsize=10)
plt.tight_layout()
final_fig = os.path.join(output_dir, f"composite_brain_maps_{measure}_global_qvalue.png")
plt.show()
plt.savefig(final_fig, dpi=300, bbox_inches='tight', facecolor='black')
plt.close(fig)
print(f"Saved composite figure for {measure}: {final_fig}")
print("\n All measure-specific brain maps saved successfully!")
# %% [markdown]
# # FDR for all predictors * measures* ROIs, FCD subtypes Storey’s q-values FDR Correction
# %% [markdown]
# # FCD subtypes group only FCDllb vs FCDlla and FCDllb vs FCDna
# %%
import os
from pathlib import Path
import logging
import pandas as pd
import numpy as np
from statsmodels.stats.multitest import multipletests
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore")
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri, r, Formula, globalenv
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
import rpy2.robjects.packages as rpackages
from rpy2.robjects.packages import importr
utils = importr('utils')
utils.chooseCRANmirror(ind=1)
# Try BiocManager if qvalue is not on CRAN
if not rpackages.isinstalled('qvalue'):
try:
if not rpackages.isinstalled('BiocManager'):
utils.install_packages('BiocManager')
biocmanager = importr('BiocManager')
biocmanager.install('qvalue')
except Exception as e:
raise RuntimeError(f"Failed to install qvalue via BiocManager: {e}")
# Now import
qvalue_pkg = importr('qvalue')
utils = importr('utils')
utils.chooseCRANmirror(ind=1)
if not rpackages.isinstalled('qvalue'):
try:
utils.install_packages('qvalue', repos='https://cloud.r-project.org', Ncpus=4)
except Exception as e:
raise RuntimeError(f"Failed to install qvalue via CRAN: {e}")
# Logging
logging.basicConfig(level=logging.INFO, format="%(asctime)s [%(levelname)s] %(message)s")
logger = logging.getLogger("betareg_global")
# R packages
utils = importr('utils')
utils.install_packages('betareg')
try:
betareg = importr('betareg')
except Exception as e:
raise RuntimeError("R package 'betareg' not available. Please install it in R.") from e
# User paths
DATA_ROOT = Path("/Volumes/groups/tohkagroup/Bonn_Epilepsy")
METADATA_XLSX = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.xlsx"
METADATA_CSV = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv"
METADATA_TSV = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv"
OUT_ROOT = DATA_ROOT / "betareg_global_results_fcdsubtypes_storey_qvalue_&_benjamini_hochberg_group_only"
OUT_ROOT.mkdir(parents=True, exist_ok=True)
# Load participant metadata
participants_df = None
if METADATA_XLSX.exists():
try:
participants_df = pd.read_excel(METADATA_XLSX, engine='openpyxl')
logger.info(f"Loaded metadata from Excel: {METADATA_XLSX}")
except Exception as e:
logger.warning(f"Failed Excel read {METADATA_XLSX}: {e}")
if participants_df is None and METADATA_CSV.exists():
try:
participants_df = pd.read_csv(METADATA_CSV, sep=';', engine='python', on_bad_lines='skip')
logger.info(f"Loaded metadata from CSV: {METADATA_CSV}")
except Exception as e:
logger.warning(f"Failed CSV read {METADATA_CSV}: {e}")
if participants_df is None and METADATA_TSV.exists():
try:
participants_df = pd.read_csv(METADATA_TSV, sep='\t', engine='python', on_bad_lines='skip')
logger.info(f"Loaded metadata from TSV: {METADATA_TSV}")
except Exception as e:
logger.warning(f"Failed TSV read {METADATA_TSV}: {e}")
if participants_df is None:
raise FileNotFoundError("Metadata not found or unreadable.")
# Clean metadata
if 'participant_id' not in participants_df.columns and 'subject_id' in participants_df.columns:
participants_df = participants_df.rename(columns={'subject_id': 'participant_id'})
# Maps
sex_map = {'F': 0, 'M': 1}
hist_map = {'IIa': 0, 'IIb': 1, 'fcdna': 2}
# Apply mappings
participants_df['histopathology'] = participants_df['histopathology'].astype(str)
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
participants_df = participants_df.dropna(subset=['participant_id', 'sex', 'histopathology', 'age_scan'])
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
logger.info(f"Loaded metadata: {len(participants_df)} participants after cleaning.")
# Settings
MEASURES = ["intensity", "thickness", "volume", "curvature", "surface_area"]
MEASURE_FILE_PATTERNS = {
"intensity": "_intensity_freesurfer_stats.csv",
"thickness": "_ThickAvg_asymmetry.csv",
"volume": "_GrayVol_asymmetry.csv",
"curvature": "_MeanCurv_asymmetry.csv",
"surface_area": "_SurfArea_asymmetry.csv",
}
GROUP_FOLDERS = ["asymmetry_fcdlla", "asymmetry_fcdllb", "asymmetry_fcdna"]
MIN_SUBJ_PER_ROI = 8
ALPHA = 0.05
CI_Z = 1.96
# Helper functions
def collect_intensity_files():
files = []
for folder in GROUP_FOLDERS:
gp = DATA_ROOT / folder
if not gp.exists():
continue
for fn in os.listdir(gp):
if fn.endswith(MEASURE_FILE_PATTERNS['intensity']):
files.append(gp / fn)
return files
def collect_nested_measure_files(measure_key):
files = []
pattern_suffix = MEASURE_FILE_PATTERNS[measure_key]
for folder in GROUP_FOLDERS:
gp = DATA_ROOT / folder
if not gp.exists():
continue
for subj in os.listdir(gp):
subj_folder = gp / subj
if not subj_folder.is_dir():
continue
stats_folder = subj_folder / ("stats" if folder=="asymmetry_hc" else "stats_cleaned")
if not stats_folder.exists():
stats_folder = subj_folder / "stats"
if not stats_folder.exists():
continue
fpath = stats_folder / f"{subj}{pattern_suffix}"
if fpath.exists():
files.append(fpath)
return files
# Main loop
all_results = []
failures = []
for measure in MEASURES:
logger.info(f"Processing measure: {measure}")
measure_out = OUT_ROOT / measure
measure_out.mkdir(parents=True, exist_ok=True)
files = collect_intensity_files() if measure=="intensity" else collect_nested_measure_files(measure)
if not files:
logger.warning(f"No files found for {measure}")
continue
rows = []
for fpath in files:
try:
df = pd.read_csv(fpath)
except:
continue
asym_col = 'Asymmetry_Index' if 'Asymmetry_Index' in df.columns else ('AsymmetryIndex' if 'AsymmetryIndex' in df.columns else None)
if asym_col is None:
continue
fname = fpath.name
participant_id = fname.replace(MEASURE_FILE_PATTERNS[measure], '') if measure=='intensity' else fname[:-len(MEASURE_FILE_PATTERNS[measure])]
df = df[['StructName', asym_col]].rename(columns={asym_col:'Asymmetry_Index'})
df['participant_id'] = participant_id
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
rows.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
if not rows:
continue
asymmetry_df = pd.concat(rows, ignore_index=True)
merged = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
merged = merged.dropna(subset=['Unsigned_Asymmetry','sex','histopathology','age_scan'])
if merged.empty:
continue
# Save distribution plot
plt.figure(figsize=(6,3))
sns.histplot(merged['Unsigned_Asymmetry'], bins=40)
plt.title(f"{measure} - Unsigned_Asymmetry distribution")
plt.tight_layout()
plt.savefig(measure_out / f"{measure}_scaled_dist.png")
plt.close()
roi_list = merged['StructName'].unique()
formula = Formula('Unsigned_Asymmetry ~ histopathology + sex + age_scan')
for roi in tqdm(roi_list, desc=f"{measure} ROIs"):
roi_df = merged[merged['StructName']==roi].copy()
# Drop rows with missing asymmetry values or missing metadata
roi_df = roi_df.dropna(subset=['Unsigned_Asymmetry','histopathology','sex','age_scan'])
roi_df = roi_df.reset_index(drop=True)
logger.info("Using predictor column: %s, unique values = %s",
'histopathology', roi_df['histopathology'].unique())
nsubj = roi_df['participant_id'].nunique()
if nsubj < MIN_SUBJ_PER_ROI or roi_df['Unsigned_Asymmetry'].nunique()<3:
failures.append({'Measure':measure,'ROI':roi,'Reason':'too_few_or_low_variability','N':nsubj})
continue
try:
with localconverter(ro.default_converter + pandas2ri.converter):
rdf = ro.conversion.py2rpy(roi_df[['Unsigned_Asymmetry','histopathology','sex','age_scan']])
globalenv['rdf'] = rdf
model = betareg.betareg(formula, data=rdf)
summary_model = r.summary(model)
# --- Residual plots ---
try:
r_resid = r.residuals(model)
with localconverter(ro.default_converter + pandas2ri.converter):
py_resid = ro.conversion.rpy2py(r_resid)
resid_series = pd.Series(py_resid).astype(float)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Histogram
axes[0].hist(resid_series, bins=30, color='skyblue', edgecolor='black')
axes[0].set_title(f"{measure} {roi} residuals histogram")
# QQ plot
import statsmodels.api as sm
sm.graphics.qqplot(resid_series, line='45', ax=axes[1])
axes[1].set_title(f"{measure} {roi} residuals QQ")
plt.tight_layout()
plt.savefig(measure_out / f"{roi}_residuals.png", dpi=200)
plt.close()
except Exception as e:
logger.warning(f"Residuals plot failed for {measure} {roi}: {e}")
coef_table = summary_model.rx2('coefficients').rx2('mean')
coef_names = list(coef_table.rownames)
for name in coef_names:
est = float(coef_table.rx(name, 'Estimate')[0])
se = float(coef_table.rx(name, 'Std. Error')[0])
z = float(coef_table.rx(name, 'z value')[0])
p = float(coef_table.rx(name, 'Pr(>|z|)')[0])
ci_low = est - CI_Z * se
ci_high = est + CI_Z * se
all_results.append({
'Measure': measure,
'ROI': roi,
'Predictor': name,
'Estimate': est,
'StdErr': se,
'z': z,
'Pval': p, # raw p-value
'CI_low': ci_low,
'CI_high': ci_high,
'Pseudo_R2': float(summary_model.rx2('pseudo.r.squared')[0]) if 'pseudo.r.squared' in summary_model.names else np.nan,
'N_subjects': nsubj
})
except Exception as e:
failures.append({'Measure':measure,'ROI':roi,'Reason':f"betareg_failed:{e}",'N':nsubj})
continue
# ---------- Results and corrections ----------
results_df = pd.DataFrame(all_results)
fails_df = pd.DataFrame(failures)
if results_df.empty:
raise RuntimeError("No model results collected.")
# --- Keep raw p-values separately for both stages ---
results_df['Pval_raw_global'] = results_df['Pval'].astype(float)
results_df['Pval_raw_predictor'] = results_df['Pval'].astype(float)
# --- Global FDR (all ROIs and measures) but only for group/histopathology ---
# Select all predictors related to histopathology
mask_interest = results_df['Predictor'].str.contains('histopathology', na=False)
pvals_all = results_df.loc[mask_interest, 'Pval_raw_global'].astype(float).values
if len(pvals_all) == 0:
raise RuntimeError("No histopathology coefficients found for global correction!")
# Benjamini-Hochberg (FDR)
fdr_adj = multipletests(pvals_all, alpha=ALPHA, method='fdr_bh')[1]
results_df.loc[mask_interest, 'FDR_p_global'] = fdr_adj
# Bonferroni correction
bonf_adj = multipletests(pvals_all, alpha=ALPHA, method='bonferroni')[1]
results_df.loc[mask_interest, 'Bonferroni_p'] = bonf_adj
# Flags
results_df.loc[mask_interest, 'significant_FDR_global'] = fdr_adj < ALPHA
results_df.loc[mask_interest, 'significant_Bonferroni'] = bonf_adj < ALPHA
# --- Storey q-values (using R's qvalue package) for global correction ---
try:
qvalue_pkg = importr('qvalue')
with localconverter(ro.default_converter + pandas2ri.converter):
r_pvals = ro.FloatVector(pvals_all)
qobj = qvalue_pkg.qvalue(r_pvals)
qvalues = np.array(qobj.rx2('qvalues'))
results_df.loc[mask_interest, 'qvalue_global'] = qvalues
results_df.loc[mask_interest, 'significant_qvalue'] = qvalues < ALPHA
logger.info("Storey q-values successfully computed for group/histopathology only (global correction).")
except Exception as e:
logger.warning(f"Storey q-value computation failed for global correction: {e}")
results_df['qvalue_global'] = np.nan
results_df['significant_qvalue'] = False
# --- Per-predictor FDR (only for 'group' or 'histopathology') ---
results_df['FDR_p_within_predictor'] = np.nan
results_df['qvalue_within_predictor'] = np.nan
results_df['significant_qvalue_within_predictor'] = False
predictors_to_correct = ['histopathology'] # <only histopathology
for (measure, predictor), grp in results_df.groupby(['Measure','Predictor']):
if 'histopathology' not in predictor:
continue # skip all other predictors like age, sex, etc.
# FDR correction (Benjamini-Hochberg)
adj_bh = multipletests(grp['Pval_raw_predictor'].values, method='fdr_bh')[1]
results_df.loc[grp.index, 'FDR_p_within_predictor'] = adj_bh
# Storey q-value
try:
with localconverter(ro.default_converter + pandas2ri.converter):
r_pvals = ro.FloatVector(grp['Pval_raw_predictor'].values)
qobj = qvalue_pkg.qvalue(r_pvals)
qvals = np.array(qobj.rx2('qvalues'))
results_df.loc[grp.index, 'qvalue_within_predictor'] = qvals
results_df.loc[grp.index, 'significant_qvalue_within_predictor'] = qvals < ALPHA
except Exception as e:
logger.warning(f"Storey q-value failed for predictor {predictor} in {measure}: {e}")
# --- Save outputs ---
results_df.to_csv(OUT_ROOT / "all_measures_betareg_results_with_global_corrections.csv", index=False)
try:
results_df.to_excel(OUT_ROOT / "all_measures_betareg_results_with_global_corrections.xlsx", index=False)
except:
logger.warning("Could not write Excel (check openpyxl)")
fails_df.to_csv(OUT_ROOT / "model_failures.csv", index=False)
logger.info(f"Saved results and failures in {OUT_ROOT}")
# %% [markdown]
# # Plot for FCD subtypes histopathology only
# %%
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from nilearn import datasets, plotting
from nibabel.freesurfer.io import read_annot
# === Directories ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/betareg_global_results_fcdsubtypes_storey_qvalue_&_benjamini_hochberg_group_only"
results_csv = os.path.join(data_root, "all_measures_betareg_results_with_global_corrections.csv")
output_dir = os.path.join(data_root, "brain_maps_fcdsubtypes_zscore_global_by_measure")
os.makedirs(output_dir, exist_ok=True)
# === Load results ===
df = pd.read_csv(results_csv)
# === Filter significant ROIs (global q < 0.05, exclude intercepts) ===
sig_df = df[
(df['qvalue_global'] < 0.05) &
(~df['Predictor'].str.lower().str.contains("intercept"))
]
if sig_df.empty:
print(" No significant ROIs found with qvalue_global < 0.05 — plotting all instead.")
sig_df = df.copy()
# === Load fsaverage5 + aparc labels ===
fsaverage = datasets.fetch_surf_fsaverage('fsaverage5')
labels_left, ctab_left, names_left = read_annot(os.path.join(data_root, "lh.aparc.annot"))
labels_right, ctab_right, names_right = read_annot(os.path.join(data_root, "rh.aparc.annot"))
names_left = [n.decode("utf-8").lower() for n in names_left]
names_right = [n.decode("utf-8").lower() for n in names_right]
# === Map ROI z-scores to surface ===
def roi_to_surface(df_pred):
# Only keep significant ROIs
roi_dict = {
roi.lower(): z
for roi, z, q in zip(df_pred['ROI'], df_pred['z'], df_pred['qvalue_global'])
if q < 0.05
}
data_left = np.zeros(len(labels_left))
data_right = np.zeros(len(labels_right))
for roi_name, z in roi_dict.items():
if roi_name in names_left:
idx = names_left.index(roi_name)
data_left[labels_left == idx] = z
if roi_name in names_right:
idx = names_right.index(roi_name)
data_right[labels_right == idx] = z
return data_left, data_right
# === Plotting function ===
def save_views(data, hemi, predictor, measure):
for view in ['lateral', 'medial']:
fname = os.path.join(output_dir, f"{measure}_{predictor}_{hemi}_{view}.png")
fig = plt.figure(facecolor='black')
ax = fig.add_subplot(111, projection='3d')
plotting.plot_surf_stat_map(
fsaverage[f'infl_{hemi}'], data,
hemi=hemi,
bg_map=fsaverage[f'sulc_{hemi}'],
cmap='bwr', colorbar=True,
bg_on_data=True, darkness=0.5, alpha=1.0, axes=ax,
vmin=-5, vmax=5, view=view
)
for child_ax in fig.get_axes():
child_ax.set_facecolor("black")
fig.savefig(fname, dpi=300, facecolor='black', bbox_inches='tight')
plt.close(fig)
# === Generate plots per measure ===
measures = ['intensity', 'thickness', 'volume', 'curvature', 'surface_area']
for measure in measures:
df_measure = sig_df[sig_df['Measure'].str.lower() == measure.lower()]
if df_measure.empty:
print(f" No significant results found for {measure}, skipping.")
continue
predictors = df_measure['Predictor'].unique()
# Step 1: generate surface images
for predictor in predictors:
df_pred = df_measure[df_measure['Predictor'] == predictor]
data_left, data_right = roi_to_surface(df_pred)
save_views(data_left, 'left', predictor, measure)
save_views(data_right, 'right', predictor, measure)
print("Saved all individual predictor brain maps to:", output_dir)
# Step 2: assemble composite for this measure
n_preds = len(predictors)
fig, axes = plt.subplots(n_preds, 4, figsize=(16, 4 * n_preds), facecolor='black')
if n_preds == 1:
axes = np.array([axes])
for i, predictor in enumerate(predictors):
for j, hemi in enumerate(['left', 'right']):
for k, view in enumerate(['lateral', 'medial']):
img_path = os.path.join(output_dir, f"{measure}_{predictor}_{hemi}_{view}.png")
if os.path.exists(img_path):
img = plt.imread(img_path)
axes[i, j*2 + k].imshow(img)
axes[i, j*2 + k].axis("off")
axes[i, j*2 + k].set_facecolor('black')
axes[i, j*2 + k].set_title(f"{predictor} ({hemi} {view})", color='yellow', fontsize=10)
plt.tight_layout()
final_fig = os.path.join(output_dir, f"composite_brain_maps_{measure}_global_qvalue.png")
plt.show()
plt.savefig(final_fig, dpi=300, bbox_inches='tight', facecolor='black')
plt.close(fig)
print(f"Saved composite figure for {measure}: {final_fig}")
print("\n All measure-specific brain maps saved successfully!")
# %% [markdown]
# # FDR for all predictors * measures* ROIs, FCDIIb vs FCDIIa Storey’s q-values FDR Correction
# %% [markdown]
# # For group or histopathology only FCDIIb vs FCDIIa
# %%
# unified_betareg_all_measures_global_corrections.py
# Single script to run beta regression across five measures, keep original file paths,
# apply global FDR + Bonferroni (FWER) across all tests,
# compute confidence intervals, and save residual plots and diagnostics.
import os
from pathlib import Path
import logging
import pandas as pd
import numpy as np
from statsmodels.stats.multitest import multipletests
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore")
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri, r, Formula, globalenv
from rpy2.robjects.packages import importr
from rpy2.robjects.conversion import localconverter
import rpy2.robjects.packages as rpackages
from rpy2.robjects.packages import importr
utils = importr('utils')
utils.chooseCRANmirror(ind=1)
# Try BiocManager if qvalue is not on CRAN
if not rpackages.isinstalled('qvalue'):
try:
if not rpackages.isinstalled('BiocManager'):
utils.install_packages('BiocManager')
biocmanager = importr('BiocManager')
biocmanager.install('qvalue')
except Exception as e:
raise RuntimeError(f"Failed to install qvalue via BiocManager: {e}")
# Now import
qvalue_pkg = importr('qvalue')
utils = importr('utils')
utils.chooseCRANmirror(ind=1)
if not rpackages.isinstalled('qvalue'):
try:
utils.install_packages('qvalue', repos='https://cloud.r-project.org', Ncpus=4)
except Exception as e:
raise RuntimeError(f"Failed to install qvalue via CRAN: {e}")
# Logging
logging.basicConfig(level=logging.INFO, format="%(asctime)s [%(levelname)s] %(message)s")
logger = logging.getLogger("betareg_global")
# R packages
utils = importr('utils')
utils.install_packages('betareg')
try:
betareg = importr('betareg')
except Exception as e:
raise RuntimeError("R package 'betareg' not available. Please install it in R.") from e
# User paths
DATA_ROOT = Path("/Volumes/groups/tohkagroup/Bonn_Epilepsy")
METADATA_XLSX = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.xlsx"
METADATA_CSV = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.csv"
METADATA_TSV = DATA_ROOT / "derivatives_fsaverage/freesurfer7.4.1/participants.aparc.tsv"
OUT_ROOT = DATA_ROOT / "betareg_global_results_fcdlla_vs_fcdllb_storey_qvalue_&_benjamini_hochberg_group_only"
OUT_ROOT.mkdir(parents=True, exist_ok=True)
# Load participant metadata
participants_df = None
if METADATA_XLSX.exists():
try:
participants_df = pd.read_excel(METADATA_XLSX, engine='openpyxl')
logger.info(f"Loaded metadata from Excel: {METADATA_XLSX}")
except Exception as e:
logger.warning(f"Failed Excel read {METADATA_XLSX}: {e}")
if participants_df is None and METADATA_CSV.exists():
try:
participants_df = pd.read_csv(METADATA_CSV, sep=';', engine='python', on_bad_lines='skip')
logger.info(f"Loaded metadata from CSV: {METADATA_CSV}")
except Exception as e:
logger.warning(f"Failed CSV read {METADATA_CSV}: {e}")
if participants_df is None and METADATA_TSV.exists():
try:
participants_df = pd.read_csv(METADATA_TSV, sep='\t', engine='python', on_bad_lines='skip')
logger.info(f"Loaded metadata from TSV: {METADATA_TSV}")
except Exception as e:
logger.warning(f"Failed TSV read {METADATA_TSV}: {e}")
if participants_df is None:
raise FileNotFoundError("Metadata not found or unreadable.")
# Clean metadata
if 'participant_id' not in participants_df.columns and 'subject_id' in participants_df.columns:
participants_df = participants_df.rename(columns={'subject_id': 'participant_id'})
participants_df['group'] = participants_df['group'].astype(str).str.lower()
participants_df['sex'] = participants_df['sex'].astype(str).str.upper()
sex_map = {'F': 0, 'M': 1}
hist_map = {'IIa': 0, 'IIb': 1}
participants_df['group'] = participants_df['histopathology'].map(hist_map)
participants_df = participants_df.dropna(subset=['participant_id', 'sex', 'group', 'age_scan'])
participants_df['participant_id'] = participants_df['participant_id'].astype(str)
logger.info(f"Loaded metadata: {len(participants_df)} participants after cleaning.")
# Settings
MEASURES = ["intensity", "thickness", "volume", "curvature", "surface_area"]
MEASURE_FILE_PATTERNS = {
"intensity": "_intensity_freesurfer_stats.csv",
"thickness": "_ThickAvg_asymmetry.csv",
"volume": "_GrayVol_asymmetry.csv",
"curvature": "_MeanCurv_asymmetry.csv",
"surface_area": "_SurfArea_asymmetry.csv",
}
GROUP_FOLDERS = ["asymmetry_fcdlla", "asymmetry_fcdllb"]
MIN_SUBJ_PER_ROI = 8
EPS = 1e-4
ALPHA = 0.05
CI_Z = 1.96
# ---------- Helper functions ----------
def collect_intensity_files():
files = []
for folder in GROUP_FOLDERS:
gp = DATA_ROOT / folder
if not gp.exists():
continue
for fn in os.listdir(gp):
if fn.endswith(MEASURE_FILE_PATTERNS['intensity']):
files.append(gp / fn)
return files
def collect_nested_measure_files(measure_key):
files = []
pattern_suffix = MEASURE_FILE_PATTERNS[measure_key]
for folder in GROUP_FOLDERS:
gp = DATA_ROOT / folder
if not gp.exists():
continue
for subj in os.listdir(gp):
subj_folder = gp / subj
if not subj_folder.is_dir():
continue
stats_folder = subj_folder / ("stats" if folder=="asymmetry_hc" else "stats_cleaned")
if not stats_folder.exists():
stats_folder = subj_folder / "stats"
if not stats_folder.exists():
continue
fpath = stats_folder / f"{subj}{pattern_suffix}"
if fpath.exists():
files.append(fpath)
return files
# ---------- Main loop ----------
all_results = []
failures = []
for measure in MEASURES:
logger.info(f"Processing measure: {measure}")
measure_out = OUT_ROOT / measure
measure_out.mkdir(parents=True, exist_ok=True)
files = collect_intensity_files() if measure=="intensity" else collect_nested_measure_files(measure)
if not files:
logger.warning(f"No files found for {measure}")
continue
rows = []
for fpath in files:
try:
df = pd.read_csv(fpath)
except:
continue
asym_col = 'Asymmetry_Index' if 'Asymmetry_Index' in df.columns else ('AsymmetryIndex' if 'AsymmetryIndex' in df.columns else None)
if asym_col is None:
continue
fname = fpath.name
participant_id = fname.replace(MEASURE_FILE_PATTERNS[measure], '') if measure=='intensity' else fname[:-len(MEASURE_FILE_PATTERNS[measure])]
df = df[['StructName', asym_col]].rename(columns={asym_col:'Asymmetry_Index'})
df['participant_id'] = participant_id
df['Unsigned_Asymmetry'] = np.abs(df['Asymmetry_Index'])
rows.append(df[['participant_id', 'StructName', 'Unsigned_Asymmetry']])
if not rows:
continue
asymmetry_df = pd.concat(rows, ignore_index=True)
merged = pd.merge(asymmetry_df, participants_df, on='participant_id', how='inner')
merged = merged.dropna(subset=['Unsigned_Asymmetry','sex','group','age_scan'])
if merged.empty:
continue
# Save distribution plot
plt.figure(figsize=(6,3))
sns.histplot(merged['Unsigned_Asymmetry'], bins=40)
plt.title(f"{measure} - Unsigned_Asymmetry distribution")
plt.tight_layout()
plt.savefig(measure_out / f"{measure}_scaled_dist.png")
plt.close()
roi_list = merged['StructName'].unique()
formula = Formula('Unsigned_Asymmetry ~ group + sex + age_scan')
for roi in tqdm(roi_list, desc=f"{measure} ROIs"):
roi_df = merged[merged['StructName'] == roi].copy()
nsubj = roi_df['participant_id'].nunique()
if nsubj < MIN_SUBJ_PER_ROI:
failures.append({'Measure': measure, 'ROI': roi, 'Reason': 'too_few_subjects', 'N': nsubj})
continue
try:
# Convert to R dataframe
with localconverter(ro.default_converter + pandas2ri.converter):
rdf = ro.conversion.py2rpy(roi_df[['Unsigned_Asymmetry', 'group', 'sex', 'age_scan']])
globalenv['rdf'] = rdf
# Fit beta regression
model = betareg.betareg(formula, data=rdf)
summary_model = r.summary(model)
# --- Residual plots ---
try:
r_resid = r.residuals(model)
with localconverter(ro.default_converter + pandas2ri.converter):
py_resid = ro.conversion.rpy2py(r_resid)
resid_series = pd.Series(py_resid).astype(float)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Histogram
axes[0].hist(resid_series, bins=30, color='skyblue', edgecolor='black')
axes[0].set_title(f"{measure} {roi} residuals histogram")
# QQ plot
import statsmodels.api as sm
sm.graphics.qqplot(resid_series, line='45', ax=axes[1])
axes[1].set_title(f"{measure} {roi} residuals QQ")
plt.tight_layout()
plt.savefig(measure_out / f"{roi}_residuals.png", dpi=200)
plt.close()
except Exception as e:
logger.warning(f"Residuals plot failed for {measure} {roi}: {e}")
# --- Coefficients table ---
coef_table = summary_model.rx2('coefficients').rx2('mean')
coef_names = list(coef_table.rownames)
for name in coef_names:
est = float(coef_table.rx(name, 'Estimate')[0])
se = float(coef_table.rx(name, 'Std. Error')[0])
z = float(coef_table.rx(name, 'z value')[0])
p = float(coef_table.rx(name, 'Pr(>|z|)')[0])
ci_low = est - CI_Z * se
ci_high = est + CI_Z * se
all_results.append({
'Measure': measure,
'ROI': roi,
'Predictor': name,
'Estimate': est,
'StdErr': se,
'z': z,
'Pval': p,
'CI_low': ci_low,
'CI_high': ci_high,
'Pseudo_R2': float(summary_model.rx2('pseudo.r.squared')[0]) if 'pseudo.r.squared' in summary_model.names else np.nan,
'N_subjects': nsubj
})
except Exception as e:
failures.append({'Measure': measure, 'ROI': roi, 'Reason': f"betareg_failed:{e}", 'N': nsubj})
continue
# ---------- Results and corrections ----------
results_df = pd.DataFrame(all_results)
fails_df = pd.DataFrame(failures)
if results_df.empty:
raise RuntimeError("No model results collected.")
# --- Keep raw p-values separately for both stages ---
results_df['Pval_raw_global'] = results_df['Pval'].astype(float)
results_df['Pval_raw_predictor'] = results_df['Pval'].astype(float)
# --- Global FDR (all ROIs and measures) but only for group/histopathology ---
mask_interest = results_df['Predictor'].isin(['group', 'histopathology'])
pvals_all = results_df.loc[mask_interest, 'Pval_raw_global'].astype(float).values
# Benjamini-Hochberg (FDR)
fdr_adj = multipletests(pvals_all, alpha=ALPHA, method='fdr_bh')[1]
results_df.loc[mask_interest, 'FDR_p_global'] = fdr_adj
# Bonferroni correction
bonf_adj = multipletests(pvals_all, alpha=ALPHA, method='bonferroni')[1]
results_df.loc[mask_interest, 'Bonferroni_p'] = bonf_adj
# Flags
results_df.loc[mask_interest, 'significant_FDR_global'] = fdr_adj < ALPHA
results_df.loc[mask_interest, 'significant_Bonferroni'] = bonf_adj < ALPHA
# --- Storey q-values (using R's qvalue package) for global correction ---
try:
qvalue_pkg = importr('qvalue')
with localconverter(ro.default_converter + pandas2ri.converter):
r_pvals = ro.FloatVector(pvals_all)
qobj = qvalue_pkg.qvalue(r_pvals)
qvalues = np.array(qobj.rx2('qvalues'))
results_df.loc[mask_interest, 'qvalue_global'] = qvalues
results_df.loc[mask_interest, 'significant_qvalue'] = qvalues < ALPHA
logger.info("Storey q-values successfully computed for group/histopathology only (global correction).")
except Exception as e:
logger.warning(f"Storey q-value computation failed for global correction: {e}")
results_df['qvalue_global'] = np.nan
results_df['significant_qvalue'] = False
logger.info("Storey q-values successfully computed and added.")
except Exception as e:
logger.warning(f"Storey q-value computation failed: {e}")
results_df['qvalue_global'] = np.nan
results_df['significant_qvalue'] = False
# --- Per-predictor FDR (only for 'group' or 'histopathology') ---
results_df['FDR_p_within_predictor'] = np.nan
results_df['qvalue_within_predictor'] = np.nan
results_df['significant_qvalue_within_predictor'] = False
predictors_to_correct = ['group', 'histopathology'] # <only group or histopathology
for (measure, predictor), grp in results_df.groupby(['Measure','Predictor']):
if predictor not in predictors_to_correct:
continue # skip all other predictors like age, sex, etc.
# FDR correction (Benjamini-Hochberg)
adj_bh = multipletests(grp['Pval_raw_predictor'].values, method='fdr_bh')[1]
results_df.loc[grp.index, 'FDR_p_within_predictor'] = adj_bh
# Storey q-value
try:
with localconverter(ro.default_converter + pandas2ri.converter):
r_pvals = ro.FloatVector(grp['Pval_raw_predictor'].values)
qobj = qvalue_pkg.qvalue(r_pvals)
qvals = np.array(qobj.rx2('qvalues'))
results_df.loc[grp.index, 'qvalue_within_predictor'] = qvals
results_df.loc[grp.index, 'significant_qvalue_within_predictor'] = qvals < ALPHA
except Exception as e:
logger.warning(f"Storey q-value failed for predictor {predictor} in {measure}: {e}")
# --- Save outputs ---
results_df.to_csv(OUT_ROOT / "all_measures_betareg_results_with_global_corrections.csv", index=False)
try:
results_df.to_excel(OUT_ROOT / "all_measures_betareg_results_with_global_corrections.xlsx", index=False)
except:
logger.warning("Could not write Excel (check openpyxl)")
fails_df.to_csv(OUT_ROOT / "model_failures.csv", index=False)
logger.info(f"Saved results and failures in {OUT_ROOT}")
# %% [markdown]
# # Plot for FCDIIb vs FCDIIa histopathology only
# %%
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from nilearn import datasets, plotting
from nibabel.freesurfer.io import read_annot
# === Directories ===
data_root = "/Volumes/groups/tohkagroup/Bonn_Epilepsy/betareg_global_results_fcdlla_vs_fcdllb_storey_qvalue_&_benjamini_hochberg_group_only"
results_csv = os.path.join(data_root, "all_measures_betareg_results_with_global_corrections.csv")
output_dir = os.path.join(data_root, "brain_maps_fcdlla_vs_fcdllb_zscore_global_by_measure")
os.makedirs(output_dir, exist_ok=True)
# === Load results ===
df = pd.read_csv(results_csv)
# === Filter significant ROIs (global q < 0.05, exclude intercepts) ===
sig_df = df[
(df['qvalue_global'] < 0.05) &
(~df['Predictor'].str.lower().str.contains("intercept"))
]
if sig_df.empty:
print(" No significant ROIs found with qvalue_global < 0.05 — plotting all instead.")
sig_df = df.copy()
# === Load fsaverage5 + aparc labels ===
fsaverage = datasets.fetch_surf_fsaverage('fsaverage5')
labels_left, ctab_left, names_left = read_annot(os.path.join(data_root, "lh.aparc.annot"))
labels_right, ctab_right, names_right = read_annot(os.path.join(data_root, "rh.aparc.annot"))
names_left = [n.decode("utf-8").lower() for n in names_left]
names_right = [n.decode("utf-8").lower() for n in names_right]
# === Map ROI z-scores to surface ===
def roi_to_surface(df_pred):
roi_dict = {roi.lower(): z for roi, z in zip(df_pred['ROI'], df_pred['z'])}
data_left = np.zeros(len(labels_left))
data_right = np.zeros(len(labels_right))
for roi_name, z in roi_dict.items():
if roi_name in names_left:
idx = names_left.index(roi_name)
data_left[labels_left == idx] = z
if roi_name in names_right:
idx = names_right.index(roi_name)
data_right[labels_right == idx] = z
return data_left, data_right
# === Plotting function ===
def save_views(data, hemi, predictor, measure):
for view in ['lateral', 'medial']:
fname = os.path.join(output_dir, f"{measure}_{predictor}_{hemi}_{view}.png")
fig = plt.figure(facecolor='black')
ax = fig.add_subplot(111, projection='3d')
plotting.plot_surf_stat_map(
fsaverage[f'infl_{hemi}'], data,
hemi=hemi,
bg_map=fsaverage[f'sulc_{hemi}'],
cmap='bwr', colorbar=True,
bg_on_data=True, darkness=0.5, alpha=1.0, axes=ax,
vmin=-5, vmax=5, view=view
)
for child_ax in fig.get_axes():
child_ax.set_facecolor("black")
fig.savefig(fname, dpi=300, facecolor='black', bbox_inches='tight')
plt.close(fig)
# === Generate plots per measure ===
measures = ['intensity', 'thickness', 'volume', 'curvature', 'surface_area']
for measure in measures:
df_measure = sig_df[sig_df['Measure'].str.lower() == measure.lower()]
if df_measure.empty:
print(f" No significant results found for {measure}, skipping.")
continue
predictors = df_measure['Predictor'].unique()
# Step 1: generate surface images
for predictor in predictors:
df_pred = df_measure[df_measure['Predictor'] == predictor]
data_left, data_right = roi_to_surface(df_pred)
save_views(data_left, 'left', predictor, measure)
save_views(data_right, 'right', predictor, measure)
print("Saved all individual predictor brain maps to:", output_dir)
# Step 2: assemble composite for this measure
n_preds = len(predictors)
fig, axes = plt.subplots(n_preds, 4, figsize=(16, 4 * n_preds), facecolor='black')
if n_preds == 1:
axes = np.array([axes])
for i, predictor in enumerate(predictors):
for j, hemi in enumerate(['left', 'right']):
for k, view in enumerate(['lateral', 'medial']):
img_path = os.path.join(output_dir, f"{measure}_{predictor}_{hemi}_{view}.png")
if os.path.exists(img_path):
img = plt.imread(img_path)
axes[i, j*2 + k].imshow(img)
axes[i, j*2 + k].axis("off")
axes[i, j*2 + k].set_facecolor('black')
axes[i, j*2 + k].set_title(f"{predictor} ({hemi} {view})", color='yellow', fontsize=10)
plt.tight_layout()
final_fig = os.path.join(output_dir, f"composite_brain_maps_{measure}_global_qvalue.png")
plt.show()
plt.savefig(final_fig, dpi=300, bbox_inches='tight', facecolor='black')
plt.close(fig)
print(f"Saved composite figure for {measure}: {final_fig}")
print("\n All measure-specific brain maps saved successfully!")
# %%
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
import os
# --- File paths ---
paths = [
"/Volumes/groups/tohkagroup/Bonn_Epilepsy/study1_fcdll_asymmetry/results/betareg_global_results_fcd_vs_hc_storey_qvalue_benjamini_hochberg_group_only/all_measures_betareg_results_with_global_corrections.xlsx",
"/Volumes/groups/tohkagroup/Bonn_Epilepsy/study1_fcdll_asymmetry/results/betareg_global_results_fcdsubtypes_storey_qvalue_benjamini_hochberg_group_only/all_measures_betareg_results_with_global_corrections.xlsx",
"/Volumes/groups/tohkagroup/Bonn_Epilepsy/study1_fcdll_asymmetry/results/betareg_global_results_fcdlla_vs_fcdllb_storey_qvalue_benjamini_hochberg_group_only/all_measures_betareg_results_with_global_corrections.xlsx"
]
# --- Corrected significance threshold ---
threshold = 3.12
# --- Load corrected z-values from the "z" column ---
z_values = []
for path in paths:
if os.path.exists(path):
df = pd.read_excel(path)
df.columns = df.columns.str.strip()
df["z"] = pd.to_numeric(df["z"].astype(str).str.replace(",", "."), errors="coerce")
z_values.extend(df["z"].dropna().values)
else:
print("File not found:", path)
z_values = np.array(z_values)
if len(z_values) == 0:
raise ValueError("No z-values loaded — check paths and column names")
# --- (based on your corrected Z distribution) ---
vmin, vmax = -4.92, 3.54
norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
cmap = plt.cm.coolwarm
cmap.set_bad(color="lightgrey") # grey non-significant
# Mask any |Z| < corrected threshold
masked = np.ma.masked_where(np.abs(z_values) < threshold, z_values)
# --- Plot colorbar only ---
fig = plt.figure(figsize=(10, 1.0))
ax = fig.add_axes([0.1, 0.4, 0.8, 0.3])
cb = plt.colorbar(
plt.cm.ScalarMappable(norm=norm, cmap=cmap),
cax=ax,
orientation="horizontal"
)
# labeling
cb.set_label("Corrected Z-values (|Z| ≥ 3.12)", fontsize=11)
cb.ax.tick_params(labelsize=10)
# Threshold markers
for x in [-threshold, threshold]:
cb.ax.axvline(x, linestyle="--", color="black", linewidth=1)
plt.show()
# %% [markdown]
# # Box plots
# %%
pip install statannotations
# %% [markdown]
# ## Unsigned Intensity Asymmetry
# %%
import pandas as pd
import os
from glob import glob
from scipy.stats import kruskal
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import matplotlib.pyplot as plt
import seaborn as sns
from statannotations.Annotator import Annotator
from matplotlib.patches import Patch
# === CONFIGURATION ===
show_as_stars = True # Set to False in the case of showing numeric p-values
base_dir = "/Volumes/groups/tohkagroup/Bonn_Epilepsy"
group_dirs = {
'hc': os.path.join(base_dir, 'asymmetry_hc'),
'fcdlla': os.path.join(base_dir, 'asymmetry_fcdlla'),
'fcdllb': os.path.join(base_dir, 'asymmetry_fcdllb'),
'fcdna': os.path.join(base_dir, 'asymmetry_fcdna'),
}
output_dir = os.path.join(base_dir, 'intensity_asymmetry_results')
os.makedirs(output_dir, exist_ok=True)
# === Load data and compute Unsigned Asymmetry ===
data_list = []
for group, path in group_dirs.items():
files = glob(os.path.join(path, '*.csv'))
for file in files:
df = pd.read_csv(file)
df['Subject_ID'] = os.path.basename(file).split('_')[0]
df['Group'] = group
# Compute unsigned asymmetry
df['Unsigned_Asymmetry'] = df['Asymmetry_Index'].abs()
data_list.append(df)
all_data = pd.concat(data_list, ignore_index=True)
# === Kruskal-Wallis ===
results = []
for struct in all_data['StructName'].unique():
subset = all_data[all_data['StructName'] == struct]
grouped = [group['Unsigned_Asymmetry'].values for name, group in subset.groupby('Group')]
if len(grouped) == len(group_dirs):
stat, p = kruskal(*grouped)
results.append({'StructName': struct, 'H-stat': stat, 'p-value': p})
results_df = pd.DataFrame(results).sort_values('p-value')
results_df['FDR'] = (results_df['p-value'] * len(results_df)).clip(upper=1.0)
# Save Kruskal-Wallis results
results_df.to_csv(os.path.join(output_dir, 'kruskal_unsigned_asymmetry.csv'), index=False)
# === Tukey HSD on top structure ===
top_struct = results_df.iloc[0]['StructName']
posthoc_data = all_data[all_data['StructName'] == top_struct]
tukey = pairwise_tukeyhsd(posthoc_data['Unsigned_Asymmetry'], posthoc_data['Group'])
# Save Tukey HSD results
tukey_df = pd.DataFrame(data=tukey._results_table.data[1:], columns=tukey._results_table.data[0])
tukey_df.to_csv(os.path.join(output_dir, f'tukey_{top_struct}_unsigned_asymmetry.csv'), index=False)
print("Top Structure:", top_struct)
print(tukey.summary())
# === Plot with annotations ===
plt.figure(figsize=(10, 6))
ax = sns.boxplot(data=posthoc_data, x='Group', y='Unsigned_Asymmetry', palette='Set2')
sns.stripplot(data=posthoc_data, x='Group', y='Unsigned_Asymmetry', color='black', alpha=0.4, jitter=True)
# Define pairs for annotation
pairs = [
("hc", "fcdlla"),
("hc", "fcdllb"),
("hc", "fcdna"),
("fcdlla", "fcdllb"),
("fcdlla", "fcdna"),
("fcdllb", "fcdna"),
]
# Build p-value dictionary from Tukey HSD summary
pval_dict = {}
for _, row in tukey_df.iterrows():
pair = tuple(sorted([row['group1'], row['group2']]))
pval_dict[pair] = float(row['p-adj'])
# Prepare p-values for annotation
annotated_pvals = [(pair, pval_dict.get(tuple(sorted(pair)), 1.0)) for pair in pairs]
pvalues = [pval for pair, pval in annotated_pvals]
# Annotate
annotator = Annotator(ax, pairs, data=posthoc_data, x='Group', y='Unsigned_Asymmetry')
annotator.set_pvalues_and_annotate(pvalues)
# Title and layout
plt.title(f"Unsigned Asymmetry Comparison - {top_struct}", fontsize=14)
plt.grid(True, linestyle='--', alpha=0.5)
# High-quality legend for significance
legend_labels = ["* p < 0.05", "** p < 0.01", "*** p < 0.001", "ns (p ≥ 0.05)"]
legend_handles = [Patch(facecolor='none', edgecolor='none', label=lbl) for lbl in legend_labels]
plt.legend(handles=legend_handles, title="Significance", loc='center left',
bbox_to_anchor=(1.02, 0.5), borderaxespad=0.5, frameon=False)
plt.tight_layout()
plt.savefig(os.path.join(output_dir, f'unsigned_asymmetry_{top_struct}.png'), dpi=300)
plt.show()
# %% [markdown]
# ## Unsigned Thickness Asymmetry
# %%
import os
import pandas as pd
import scipy.stats as stats
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from statannotations.Annotator import Annotator
# Define the output directory
output_dir = '/Volumes/groups/tohkagroup/Bonn_Epilepsy/thick_asymetry_result'
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# Group directories
base_dirs = {
'hc': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_hc',
'fcdlla': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdlla',
'fcdllb': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdllb',
'fcdna': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna',
}
all_data = []
print("Searching for valid asymmetry files...")
for group, base_path in base_dirs.items():
for subj_name in os.listdir(base_path):
if not subj_name.startswith('sub-'):
continue
subject_path = os.path.join(base_path, subj_name)
stats_folder = 'stats' if group == 'hc' else 'stats_cleaned'
stats_path = os.path.join(subject_path, stats_folder)
if not os.path.isdir(stats_path):
continue
# Find correct *_ThickAvg_asymmetry.csv file
for file in os.listdir(stats_path):
if file.endswith('_ThickAvg_asymmetry.csv') and file.startswith(subj_name):
full_path = os.path.join(stats_path, file)
try:
df = pd.read_csv(full_path)
df['Subject'] = subj_name
df['Group'] = group
all_data.append(df[['Subject', 'Group', 'StructName', 'Asymmetry_Index']])
except Exception as e:
print(f"Error reading {full_path}: {e}")
break # Use only the first matching file
# Combine into a single dataframe
asym_df = pd.concat(all_data, ignore_index=True)
# Check for missing values in the Asymmetry_Index column
asym_df = asym_df.dropna(subset=['Asymmetry_Index'])
# Convert to absolute asymmetry (unsigned)
asym_df['Asymmetry_Index'] = asym_df['Asymmetry_Index'].abs()
print(f"\nLoaded data from {asym_df['Subject'].nunique()} unique subjects.")
print(f"Total brain structures: {asym_df['StructName'].nunique()}.\n")
# Function to calculate Cohen's d (effect size)
def cohen_d(group1, group2):
pooled_std = (((len(group1) - 1) * group1.std()**2 + (len(group2) - 1) * group2.std()**2) /
(len(group1) + len(group2) - 2))**0.5
return (group1.mean() - group2.mean()) / pooled_std
# === Kruskal-Wallis Test (Non-parametric ANOVA) ===
kruskal_results = []
for struct in asym_df['StructName'].unique():
subset = asym_df[asym_df['StructName'] == struct]
grouped = [group['Asymmetry_Index'].values for name, group in subset.groupby('Group')]
if len(grouped) == 4: # All groups present
stat, p = stats.kruskal(*grouped)
kruskal_results.append({'StructName': struct, 'H-stat': stat, 'p-value': p})
kruskal_df = pd.DataFrame(kruskal_results)
kruskal_df['FDR'] = multipletests(kruskal_df['p-value'], method='fdr_bh')[1]
kruskal_df = kruskal_df.sort_values('p-value')
# Save Kruskal-Wallis results if needed
kruskal_csv_path = os.path.join(output_dir, 'kruskal_results.csv')
kruskal_df.to_csv(kruskal_csv_path, index=False)
print(f"Kruskal-Wallis results saved to: {kruskal_csv_path}")
# === Pairwise t-tests: HC vs each FCD group ===
t_test_results = []
for struct in asym_df['StructName'].unique():
struct_data = asym_df[asym_df['StructName'] == struct]
hc_values = struct_data[struct_data['Group'] == 'hc']['Asymmetry_Index']
for fcd_group in ['fcdlla', 'fcdllb', 'fcdna']:
fcd_values = struct_data[struct_data['Group'] == fcd_group]['Asymmetry_Index']
if len(hc_values) > 1 and len(fcd_values) > 1:
t_stat, p_val = stats.ttest_ind(hc_values, fcd_values, equal_var=False) # Welch's t-test
effect_size = cohen_d(hc_values, fcd_values) # Cohen's d for effect size
t_test_results.append({
'StructName': struct,
'Comparison': f'hc vs {fcd_group}',
't-stat': t_stat,
'p-value': p_val,
'cohen_d': effect_size
})
t_test_df = pd.DataFrame(t_test_results)
t_test_df['FDR'] = multipletests(t_test_df['p-value'], method='fdr_bh')[1]
t_test_df = t_test_df.sort_values('p-value')
# Save t-test results if needed
t_test_csv_path = os.path.join(output_dir, 't_test_results.csv')
t_test_df.to_csv(t_test_csv_path, index=False)
print(f"T-test results saved to: {t_test_csv_path}")
# === Optional: Tukey's HSD for top structure ===
top_struct = kruskal_df.iloc[0]['StructName']
posthoc_data = asym_df[asym_df['StructName'] == top_struct]
tukey = pairwise_tukeyhsd(posthoc_data['Asymmetry_Index'], posthoc_data['Group'])
print(f"Top structure by Kruskal-Wallis: {top_struct}")
print(tukey)
# === Optional: Boxplot visualization ===
plt.figure(figsize=(10, 6))
ax = sns.boxplot(data=posthoc_data, x='Group', y='Asymmetry_Index', palette='Set2')
sns.stripplot(data=posthoc_data, x='Group', y='Asymmetry_Index', color='black', alpha=0.3, jitter=True)
# Add statistical annotations
pairs = [
("hc", "fcdlla"),
("hc", "fcdllb"),
("hc", "fcdna"),
("fcdlla", "fcdllb"),
("fcdlla", "fcdna"),
("fcdllb", "fcdna")
]
# Build p-value dictionary from Tukey HSD summary
tukey_df = pd.DataFrame(data=tukey._results_table.data[1:], columns=tukey._results_table.data[0])
pval_dict = {}
for _, row in tukey_df.iterrows():
pair = tuple(sorted([row['group1'], row['group2']]))
pval_dict[pair] = float(row['p-adj'])
# Annotate p-values with significance stars
annotated_pvals = [(pair, pval_dict.get(tuple(sorted(pair)), 1.0)) for pair in pairs]
pvalues = [pval for pair, pval in annotated_pvals]
# Annotator object for stars
annotator = Annotator(ax, pairs, data=posthoc_data, x='Group', y='Asymmetry_Index')
annotator.set_pvalues_and_annotate(pvalues) # Add stars
# Title and grid
plt.title(f"Gray Matter Thickness Asymmetry Index Comparison of {top_struct}")
plt.grid(True)
plt.tight_layout()
# Save the plot as PNG
plot_filename = f"{top_struct}_asymmetry_plot.png"
plt.savefig(os.path.join(output_dir, plot_filename))
# High-quality legend for significance levels
legend_labels = [
"* p < 0.05",
"** p < 0.01",
"*** p < 0.001",
"ns (p ≥ 0.05)"
]
legend_handles = [Patch(facecolor='none', edgecolor='none', label=lbl) for lbl in legend_labels]
plt.legend(handles=legend_handles, title="Significance", loc='center left',
bbox_to_anchor=(1.02, 0.5), borderaxespad=0.5, frameon=False)
# Display the plot
plt.show()
plt.close()
print(f"Plot saved to: {os.path.join(output_dir, plot_filename)}")
# === Display the content of the CSV files ===
print("\nContent of the Kruskal-Wallis results:")
print(kruskal_df.head())
print("\nContent of the t-test results:")
print(t_test_df.head())
print("\nDone.")
# %% [markdown]
# ## Unsigned Gray Volume Asymmetry
# %%
import os
import pandas as pd
import scipy.stats as stats
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from statannotations.Annotator import Annotator
# Define the output directory
output_dir = '/Volumes/groups/tohkagroup/Bonn_Epilepsy/gvol_asymetry_result'
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# Group directories
base_dirs = {
'hc': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_hc',
'fcdlla': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdlla',
'fcdllb': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdllb',
'fcdna': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna',
}
all_data = []
print("Searching for valid asymmetry files...")
for group, base_path in base_dirs.items():
for subj_name in os.listdir(base_path):
if not subj_name.startswith('sub-'):
continue
subject_path = os.path.join(base_path, subj_name)
stats_folder = 'stats' if group == 'hc' else 'stats_cleaned'
stats_path = os.path.join(subject_path, stats_folder)
if not os.path.isdir(stats_path):
continue
# Find correct *_GrayVol_asymmetry.csv file
for file in os.listdir(stats_path):
if file.endswith('_GrayVol_asymmetry.csv') and file.startswith(subj_name):
full_path = os.path.join(stats_path, file)
try:
df = pd.read_csv(full_path)
df['Subject'] = subj_name
df['Group'] = group
all_data.append(df[['Subject', 'Group', 'StructName', 'Asymmetry_Index']])
except Exception as e:
print(f"Error reading {full_path}: {e}")
break # Use only the first matching file
# Combine into a single dataframe
asym_df = pd.concat(all_data, ignore_index=True)
# Check for missing values in the Asymmetry_Index column
asym_df = asym_df.dropna(subset=['Asymmetry_Index'])
# Convert to absolute asymmetry (unsigned)
asym_df['Asymmetry_Index'] = asym_df['Asymmetry_Index'].abs()
print(f"\nLoaded data from {asym_df['Subject'].nunique()} unique subjects.")
print(f"Total brain structures: {asym_df['StructName'].nunique()}.\n")
# Function to calculate Cohen's d (effect size)
def cohen_d(group1, group2):
pooled_std = (((len(group1) - 1) * group1.std()**2 + (len(group2) - 1) * group2.std()**2) /
(len(group1) + len(group2) - 2))**0.5
return (group1.mean() - group2.mean()) / pooled_std
# === Kruskal-Wallis Test (Non-parametric ANOVA) ===
kruskal_results = []
for struct in asym_df['StructName'].unique():
subset = asym_df[asym_df['StructName'] == struct]
grouped = [group['Asymmetry_Index'].values for name, group in subset.groupby('Group')]
if len(grouped) == 4: # All groups present
stat, p = stats.kruskal(*grouped)
kruskal_results.append({'StructName': struct, 'H-stat': stat, 'p-value': p})
kruskal_df = pd.DataFrame(kruskal_results)
kruskal_df['FDR'] = multipletests(kruskal_df['p-value'], method='fdr_bh')[1]
kruskal_df = kruskal_df.sort_values('p-value')
# Save Kruskal-Wallis results if needed
kruskal_csv_path = os.path.join(output_dir, 'kruskal_results.csv')
kruskal_df.to_csv(kruskal_csv_path, index=False)
print(f"Kruskal-Wallis results saved to: {kruskal_csv_path}")
# === Pairwise t-tests: HC vs each FCD group ===
t_test_results = []
for struct in asym_df['StructName'].unique():
struct_data = asym_df[asym_df['StructName'] == struct]
hc_values = struct_data[struct_data['Group'] == 'hc']['Asymmetry_Index']
for fcd_group in ['fcdlla', 'fcdllb', 'fcdna']:
fcd_values = struct_data[struct_data['Group'] == fcd_group]['Asymmetry_Index']
if len(hc_values) > 1 and len(fcd_values) > 1:
t_stat, p_val = stats.ttest_ind(hc_values, fcd_values, equal_var=False) # Welch's t-test
effect_size = cohen_d(hc_values, fcd_values) # Cohen's d for effect size
t_test_results.append({
'StructName': struct,
'Comparison': f'hc vs {fcd_group}',
't-stat': t_stat,
'p-value': p_val,
'cohen_d': effect_size
})
t_test_df = pd.DataFrame(t_test_results)
t_test_df['FDR'] = multipletests(t_test_df['p-value'], method='fdr_bh')[1]
t_test_df = t_test_df.sort_values('p-value')
# Save t-test results if needed
t_test_csv_path = os.path.join(output_dir, 't_test_results.csv')
t_test_df.to_csv(t_test_csv_path, index=False)
print(f"T-test results saved to: {t_test_csv_path}")
# === Optional: Tukey's HSD for top structure ===
top_struct = kruskal_df.iloc[0]['StructName']
posthoc_data = asym_df[asym_df['StructName'] == top_struct]
tukey = pairwise_tukeyhsd(posthoc_data['Asymmetry_Index'], posthoc_data['Group'])
print(f"Top structure by Kruskal-Wallis: {top_struct}")
print(tukey)
# === Optional: Boxplot visualization ===
plt.figure(figsize=(10, 6))
ax = sns.boxplot(data=posthoc_data, x='Group', y='Asymmetry_Index', palette='Set2')
sns.stripplot(data=posthoc_data, x='Group', y='Asymmetry_Index', color='black', alpha=0.3, jitter=True)
# Add statistical annotations
pairs = [
("hc", "fcdlla"),
("hc", "fcdllb"),
("hc", "fcdna"),
("fcdlla", "fcdllb"),
("fcdlla", "fcdna"),
("fcdllb", "fcdna")
]
# Build p-value dictionary from Tukey HSD summary
tukey_df = pd.DataFrame(data=tukey._results_table.data[1:], columns=tukey._results_table.data[0])
pval_dict = {}
for _, row in tukey_df.iterrows():
pair = tuple(sorted([row['group1'], row['group2']]))
pval_dict[pair] = float(row['p-adj'])
# Annotate p-values with significance stars
annotated_pvals = [(pair, pval_dict.get(tuple(sorted(pair)), 1.0)) for pair in pairs]
pvalues = [pval for pair, pval in annotated_pvals]
# Annotator object for stars
annotator = Annotator(ax, pairs, data=posthoc_data, x='Group', y='Asymmetry_Index')
annotator.set_pvalues_and_annotate(pvalues) # Add stars
# Title and grid
plt.title(f"Gray Matter Volume Asymmetry Index Comparison of {top_struct}")
plt.grid(True)
plt.tight_layout()
# Save the plot as PNG
plot_filename = f"{top_struct}_asymmetry_plot.png"
plt.savefig(os.path.join(output_dir, plot_filename))
# High-quality legend for significance levels
legend_labels = [
"* p < 0.05",
"** p < 0.01",
"*** p < 0.001",
"ns (p ≥ 0.05)"
]
legend_handles = [Patch(facecolor='none', edgecolor='none', label=lbl) for lbl in legend_labels]
plt.legend(handles=legend_handles, title="Significance", loc='center left',
bbox_to_anchor=(1.02, 0.5), borderaxespad=0.5, frameon=False)
# Display the plot
plt.show()
plt.close()
print(f"Plot saved to: {os.path.join(output_dir, plot_filename)}")
# === Display the content of the CSV files ===
print("\nContent of the Kruskal-Wallis results:")
print(kruskal_df.head())
print("\nContent of the t-test results:")
print(t_test_df.head())
print("\nDone.")
# %% [markdown]
# ## Unsigned Curvature Asymmetry
# %%
import os
import pandas as pd
import scipy.stats as stats
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from statannotations.Annotator import Annotator
# Define the output directory
output_dir = '/Volumes/groups/tohkagroup/Bonn_Epilepsy/curv_asymetry_result'
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# Group directories
base_dirs = {
'hc': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_hc',
'fcdlla': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdlla',
'fcdllb': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdllb',
'fcdna': '/Volumes/groups/tohkagroup/Bonn_Epilepsy/asymmetry_fcdna',
}
all_data = []
print("Searching for valid asymmetry files...")
for group, base_path in base_dirs.items():
for subj_name in os.listdir(base_path):
if not subj_name.startswith('sub-'):
continue
subject_path = os.path.join(base_path, subj_name)
stats_folder = 'stats' if group == 'hc' else 'stats_cleaned'
stats_path = os.path.join(subject_path, stats_folder)
if not os.path.isdir(stats_path):
continue
# Find correct *_MeanCurv_asymmetry.csv file
for file in os.listdir(stats_path):
if file.endswith('_MeanCurv_asymmetry.csv') and file.startswith(subj_name):
full_path = os.path.join(stats_path, file)
try:
df = pd.read_csv(full_path)
df['Subject'] = subj_name
df['Group'] = group
all_data.append(df[['Subject', 'Group', 'StructName', 'Asymmetry_Index']])
except Exception as e:
print(f"Error reading {full_path}: {e}")
break # Use only the first matching file
# Combine into a single dataframe
asym_df = pd.concat(all_data, ignore_index=True)
# Check for missing values in the Asymmetry_Index column
asym_df = asym_df.dropna(subset=['Asymmetry_Index'])
# Convert to absolute asymmetry (unsigned)
asym_df['Asymmetry_Index'] = asym_df['Asymmetry_Index'].abs()
print(f"\nLoaded data from {asym_df['Subject'].nunique()} unique subjects.")
print(f"Total brain structures: {asym_df['StructName'].nunique()}.\n")
# Function to calculate Cohen's d (effect size)
def cohen_d(group1, group2):
pooled_std = (((len(group1) - 1) * group1.std()**2 + (len(group2) - 1) * group2.std()**2) /
(len(group1) + len(group2) - 2))**0.5
return (group1.mean() - group2.mean()) / pooled_std
# === Kruskal-Wallis Test (Non-parametric ANOVA) ===
kruskal_results = []
for struct in asym_df['StructName'].unique():
subset = asym_df[asym_df['StructName'] == struct]
grouped = [group['Asymmetry_Index'].values for name, group in subset.groupby('Group')]
if len(grouped) == 4: # All groups present
stat, p = stats.kruskal(*grouped)
kruskal_results.append({'StructName': struct, 'H-stat': stat, 'p-value': p})
kruskal_df = pd.DataFrame(kruskal_results)
kruskal_df['FDR'] = multipletests(kruskal_df['p-value'], method='fdr_bh')[1]
kruskal_df = kruskal_df.sort_values('p-value')
# Save Kruskal-Wallis results if needed
kruskal_csv_path = os.path.join(output_dir, 'kruskal_results.csv')
kruskal_df.to_csv(kruskal_csv_path, index=False)
print(f"Kruskal-Wallis results saved to: {kruskal_csv_path}")
# === Pairwise t-tests: HC vs each FCD group ===
t_test_results = []
for struct in asym_df['StructName'].unique():
struct_data = asym_df[asym_df['StructName'] == struct]
hc_values = struct_data[struct_data['Group'] == 'hc']['Asymmetry_Index']
for fcd_group in ['fcdlla', 'fcdllb', 'fcdna']:
fcd_values = struct_data[struct_data['Group'] == fcd_group]['Asymmetry_Index']
if len(hc_values) > 1 and len(fcd_values) > 1:
t_stat, p_val = stats.ttest_ind(hc_values, fcd_values, equal_var=False) # Welch's t-test
effect_size = cohen_d(hc_values, fcd_values) # Cohen's d for effect size
t_test_results.append({
'StructName': struct,
'Comparison': f'hc vs {fcd_group}',
't-stat': t_stat,
'p-value': p_val,
'cohen_d': effect_size
})
t_test_df = pd.DataFrame(t_test_results)
t_test_df['FDR'] = multipletests(t_test_df['p-value'], method='fdr_bh')[1]
t_test_df = t_test_df.sort_values('p-value')
# Save t-test results if needed
t_test_csv_path = os.path.join(output_dir, 't_test_results.csv')
t_test_df.to_csv(t_test_csv_path, index=False)
print(f"T-test results saved to: {t_test_csv_path}")
# === Optional: Tukey's HSD for top structure ===
top_struct = kruskal_df.iloc[0]['StructName']
posthoc_data = asym_df[asym_df['StructName'] == top_struct]
tukey = pairwise_tukeyhsd(posthoc_data['Asymmetry_Index'], posthoc_data['Group'])
print(f"Top structure by Kruskal-Wallis: {top_struct}")
print(tukey)
# === Optional: Boxplot visualization ===
plt.figure(figsize=(10, 6))
ax = sns.boxplot(data=posthoc_data, x='Group', y='Asymmetry_Index', palette='Set2')
sns.stripplot(data=posthoc_data, x='Group', y='Asymmetry

Bonn_visualization.ipynb at commit bcc3292, under MIT · at the source

Overview

Authors: Faezeh Heidari¹, Mastaneh Torkamani‐Azar², Jussi Tohka², Reetta Kälviäinen^1,3

ORCID iDs: Faezeh Heidari, Mastaneh Torkamani‐Azar, Jussi Tohka, Reetta Kälviäinen

Institutions: University of Eastern Finland (Finland); Kuopio University Hospital (Finland)

Journal: Epilepsia open, article 10.1002/epi4.70336

Dates: received 12 March 2026; accepted 27 July 2026; published online 5 September 2026; in print September 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1002/epi4.70336 · PMID 42700155 · PMCID PMC13545978 · OpenAlex W7208828703

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: structural MRI / diffusion (modality), human (organism), epilepsy (population), clinical / translational (subfield)

Methods: Connectivity, Statistics, fMRI & imaging, Preprocessing

Keywords: cortex, epilepsy, focal cortical dysplasia, interhemispheric asymmetry

Topic: Epilepsy research and treatment (Psychiatry and Mental health, Medicine), according to OpenAlex

Funding: Research Council of Finland (358944); Jane ja Aatos Erkon Säätiö; European Union’s Horizon 2020 (101034307)

Citations: not cited yet (Europe PMC); 44 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.

Repository

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

faezeheidari/FCD_Asymmetry

License: MIT

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: bcc32922fc78cf2afd6990a0e0a54a7b9520eab8, 6 March 2026

Languages: Jupyter (1)

Size: 7 files, 1 script

Software Heritage: not archived

Found in: “CODE AVAILABILITY STATEMENT”

Holds: README, license file, 1 notebook

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

Tools: Matplotlib (1 file), NiBabel (1 file), Nilearn (1 file), NumPy (1 file), pandas (1 file), rpy2 (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statannotations (1 file), statsmodels (1 file)

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

26 September 2026: the link answers

3 files

Bonn_visualization.ipynb, Jupyter, 5,297 lines, 2 matches
LICENSE, License, 21 lines
README.md, Text, 69 lines

At the source: github.com/faezeheidari/FCD_Asymmetry

Code availability statement

The paper has a code 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: faezeheidari/FCD_Asymmetry

Read it in the paper: doi.org/10.1002/epi4.70336.

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;
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

Datasets cited

openneuro:ds004199, at OpenNeuro; found in “DATA AVAILABILITY STATEMENT”

Data availability statement

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

it points to a dataset: OpenNeuro ds004199
it says that the data are available on request

Read it in the paper: doi.org/10.1002/epi4.70336.

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, pages, dates, 4 authors, 4 keywords, 3 funders, 41 references.

Cite

This paper

Heidari, F., Torkamani‐Azar, M., Tohka, J., & Kälviäinen, R. (2026). Quantifying extra-lesional interhemispheric cortical asymmetry in focal cortical dysplasia type II. Epilepsia open, 10.1002/epi4.70336. https://doi.org/10.1002/epi4.70336

BibTeX

@article{heidari2026quantifying, author = {Heidari, Faezeh and Torkamani‐Azar, Mastaneh and Tohka, Jussi and Kälviäinen, Reetta}, title = {{Quantifying extra-lesional interhemispheric cortical asymmetry in focal cortical dysplasia type II}}, journal = {Epilepsia open}, year = {2026}, month = sep, pages = {10.1002/epi4.70336}, publisher = {Wiley}, issn = {2470-9239}, doi = {10.1002/epi4.70336}, url = {https://doi.org/10.1002/epi4.70336}, pmid = {42700155}, pmcid = {PMC13545978} }

RIS

TY - JOUR AU - Heidari, Faezeh AU - Torkamani‐Azar, Mastaneh AU - Tohka, Jussi AU - Kälviäinen, Reetta TI - Quantifying extra-lesional interhemispheric cortical asymmetry in focal cortical dysplasia type II T2 - Epilepsia open J2 - Epilepsia Open PY - 2026 DA - 2026/09/05 SP - 10.1002/epi4.70336 SN - 2470-9239 PB - Wiley DO - 10.1002/epi4.70336 UR - https://doi.org/10.1002/epi4.70336 LA - en ER -

CSL-JSON

{ "id": "10.1002/epi4.70336", "type": "article-journal", "title": "Quantifying extra-lesional interhemispheric cortical asymmetry in focal cortical dysplasia type II", "container-title": "Epilepsia open", "author": [ { "family": "Heidari", "given": "Faezeh" }, { "family": "Torkamani‐Azar", "given": "Mastaneh" }, { "family": "Tohka", "given": "Jussi" }, { "family": "Kälviäinen", "given": "Reetta" } ], "container-title-short": "Epilepsia Open", "page": "10.1002/epi4.70336", "DOI": "10.1002/epi4.70336", "PMID": "42700155", "PMCID": "PMC13545978", "ISSN": "2470-9239", "publisher": "Wiley", "URL": "https://doi.org/10.1002/epi4.70336", "language": "en", "issued": { "date-parts": [ [ 2026, 9, 5 ] ] } }

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.1007/s00234-026-04103-8 [code]: Enhanced detection of subtle cortical abnormalities in focal epilepsy using 7 T MRI surface-based models and graph neural networks.
Journal: Neuroradiology
In common: statannotations, Nilearn, statsmodels, 7 other tools, epilepsy, structural MRI / diffusion, 7 references
[2] doi:10.1016/j.isci.2026.117446 [code]: Volume-inflation registration (INFREG) for morphometric analysis of human focal cortical dysplasia type II.
Journal: iScience
In common: NiBabel, pandas, SciPy, 2 other tools, OpenNeuro ds004199, epilepsy, structural MRI / diffusion, 5 references
[3] doi:10.64898/2026.08.18.26360725 [code]: Temporal pole blurring in hippocampal sclerosis reflects seizure-disrupted myelination
Journal: medRxiv (preprint)
In common: statannotations, Nilearn, statsmodels, 7 other tools, epilepsy, structural MRI / diffusion, 3 references
[4] doi:10.64898/2026.04.02.26349812 [code]: 9.4 Tesla MRI in focal epilepsy patients with high-resolution surface-based profiling of focal cortical dysplasias
Journal: medRxiv (preprint)
In common: Nilearn, NiBabel, seaborn, 4 other tools, epilepsy, structural MRI / diffusion, clinical / translational, 5 references
[5] doi:10.1186/s40708-026-00299-w: Pipeline evaluation of a state-of-the-art AI algorithm for detection of focal cortical dysplasia: insights into potential failure sources.
Journal: Brain informatics
In common: OpenNeuro ds004199, epilepsy, structural MRI / diffusion, 6 references
[6] doi:10.1126/sciadv.adu9309 [code]: Variations of global brain asymmetry are associated with aging and related diseases.
Journal: Science advances
In common: Nilearn, NiBabel, scikit-learn, 3 other tools, 5 references
[7] doi:10.1371/journal.pbio.3003856 [code]: Aging and metabolism contribute separately to brain-body health.
Journal: PLoS biology
In common: rpy2, Nilearn, statsmodels, 7 other tools, structural MRI / diffusion, clinical / translational, 1 reference
[8] doi:10.21203/rs.3.rs-9914920/v1 [code]: Prediction of cognitive performance by demographics, sleep, and brain morphometry: machine learning findings from ENIGMA-Sleep Working Group
Journal: Research Square (preprint)
In common: statannotations, Nilearn, statsmodels, 7 other tools, structural MRI / diffusion
[9] doi:10.1038/s41598-026-51531-w [code]: Multimodal age-dependent diffusion-MRI analysis of the neocortex in a rat model of cortical dysplasia.
Journal: Scientific reports
In common: NiBabel, scikit-learn, pandas, 3 other tools, structural MRI / diffusion, 4 references
[10] doi:10.1162/imag.a.1256 [code]: Gamer in the scanner: Event-related analysis of fMRI activity during retro videogame play guided by automated annotations of game content.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: statannotations, Nilearn, statsmodels, 7 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.