Morphometrics Analysis - Macaque Corpus Callosum G-Ratio Analysis

This notebook loads per-axon morphometrics computed by AxonDeepSeg across 8 regions of the macaque corpus callosum, applies quality filters, and generates publication-quality plots of:

Axon diameter vs Fiber diameter with linear regression
G-ratio vs Fiber diameter with linear fit

Setup¶

import pandas as pd
import numpy as np
from scipy import stats as scipy_stats
from pathlib import Path
from IPython.display import display
from PIL import Image
import plotly.graph_objects as go
from plotly.subplots import make_subplots

DATA_DIR = Path("data/macaque/data")

SAMPLE_SLICES = {
    "sample1": ["01", "03", "07"],
    "sample2": ["01", "03", "05"],
    "sample3": ["01", "03", "05"],
    "sample4": ["01", "05", "07"],
    "sample5": ["01", "03", "05"],
    "sample6": ["01", "03", "05"],
    "sample7": ["01", "03", "05"],
    "sample8": ["01", "03", "05"],
}

REGION_LABELS = {
    "sample1": "Genu (1)",
    "sample2": "(2)",
    "sample3": "(3)",
    "sample4": "(4)",
    "sample5": "(5)",
    "sample6": "(6)",
    "sample7": "(7)",
    "sample8": "Splenium (8)",
}

COLUMN_RENAME = {
    "x0 (px)": "x0", "y0 (px)": "y0",
    "axon_area (um^2)": "axon_area", "axon_perimeter (um)": "axon_perimeter",
    "myelin_area (um^2)": "myelin_area", "axon_diam (um)": "axon_diam",
    "myelin_thickness (um)": "myelin_thickness",
    "axonmyelin_area (um^2)": "axonmyelin_area",
    "axonmyelin_perimeter (um)": "axonmyelin_perimeter",
}

# Reference values from Stikov et al. 2015
# P1: Stikov et al., NeuroImage 118, 397-405 (2015) — MRI g-ratio and aggregate histology g-ratio
# P2: Stikov et al., Data in Brief 4, 368-373 (2015) — per-axon morphometrics from slice 01 of each region
PAPER_REF = {
    "sample1": {"mri_g": "0.63\u00b10.05", "agg_g_p1": "0.69\u00b10.05", "agg_g_p2": "0.70", "mean_g_p2": "0.65"},
    "sample2": {"mri_g": "0.68\u00b10.04", "agg_g_p1": "0.68\u00b10.04", "agg_g_p2": "0.72", "mean_g_p2": "0.66"},
    "sample3": {"mri_g": "0.68\u00b10.04", "agg_g_p1": "0.64\u00b10.06", "agg_g_p2": "0.67", "mean_g_p2": "0.62"},
    "sample4": {"mri_g": "0.76\u00b10.03", "agg_g_p1": "0.69\u00b10.04", "agg_g_p2": "0.75", "mean_g_p2": "0.68"},
    "sample5": {"mri_g": "0.79\u00b10.04", "agg_g_p1": "0.72\u00b10.03", "agg_g_p2": "0.72", "mean_g_p2": "0.65"},
    "sample6": {"mri_g": "0.72\u00b10.06", "agg_g_p1": "0.64\u00b10.01", "agg_g_p2": "0.74", "mean_g_p2": "0.66"},
    "sample7": {"mri_g": "0.63\u00b10.04", "agg_g_p1": "0.60\u00b10.04", "agg_g_p2": "0.67", "mean_g_p2": "0.60"},
    "sample8": {"mri_g": "0.78\u00b10.03", "agg_g_p1": "0.74\u00b10.06", "agg_g_p2": "0.75", "mean_g_p2": "0.65"},
}

Load and filter morphometrics¶

Each morphometrics file (*_axon_morphometrics.xlsx) was generated by axondeepseg_morphometrics using manually corrected segmentation masks and a pixel size of 0.009144 μm.

The following filters are applied to each slice:

Border-touching axons are removed (their measurements are unreliable due to truncation).
Unmyelinated axons are removed (myelin area ≤ 0, g-ratio ≥ 1.0, or myelin thickness ≤ 0).
Invalid detections are removed (axon area ≤ 0 or axon diameter ≤ 0).
Rows with NaN values in key columns are dropped.

def load_morphometrics(sample, slice_id):
    """Load a single morphometrics xlsx file."""
    path = DATA_DIR / sample / f"{slice_id}_axon_morphometrics.xlsx"
    df = pd.read_excel(path)
    df = df.rename(columns=COLUMN_RENAME)
    df["sample"] = sample
    df["slice_id"] = slice_id
    return df


def filter_morphometrics(df):
    """Filter out invalid axons. Returns filtered DataFrame."""
    if "image_border_touching" in df.columns:
        df = df[df["image_border_touching"] != True]
    df = df[df["myelin_area"] > 0]
    df = df[df["myelin_thickness"] > 0]
    df = df[df["gratio"] < 1.0]
    df = df[df["axon_area"] > 0]
    df = df[df["axon_diam"] > 0]
    df = df[df["gratio"] > 0]
    df = df.dropna(subset=["gratio", "axon_diam", "myelin_thickness", "myelin_area"])
    return df

# Load, filter, and combine morphometrics for each sample
# Also compute aggregate g-ratio per image from binary masks:
#   AVF = axon_pixels / total_pixels
#   MVF = myelin_pixels / total_pixels
#   g_agg = sqrt(1 - MVF/FVF) = sqrt(AVF / (AVF + MVF))
all_data = {}
agg_gratio = {}  # sample -> list of per-image aggregate g-ratios
filter_summary = []

for sample in sorted(SAMPLE_SLICES.keys()):
    dfs = []
    sample_agg = []
    for slice_id in SAMPLE_SLICES[sample]:
        df = load_morphometrics(sample, slice_id)
        n_before = len(df)
        df = filter_morphometrics(df)
        n_after = len(df)
        filter_summary.append({
            "Region": REGION_LABELS[sample],
            "Slice": slice_id,
            "Before": n_before,
            "After": n_after,
            "Removed": n_before - n_after,
        })
        dfs.append(df)

        # Aggregate g-ratio from binary masks
        axon_mask = np.array(Image.open(DATA_DIR / sample / f"{slice_id}_seg-axon.png"))
        myelin_mask = np.array(Image.open(DATA_DIR / sample / f"{slice_id}_seg-myelin.png"))
        avf = np.sum(axon_mask > 0)
        mvf = np.sum(myelin_mask > 0)
        if avf + mvf > 0:
            g_agg = np.sqrt(avf / (avf + mvf))
        else:
            g_agg = np.nan
        sample_agg.append(g_agg)

    combined = pd.concat(dfs, ignore_index=True)
    combined["fiber_diam"] = combined["axon_diam"] + 2 * combined["myelin_thickness"]
    all_data[sample] = combined
    agg_gratio[sample] = sample_agg

display(pd.DataFrame(filter_summary))

Summary statistics¶

The aggregate g-ratio is computed per image from the binary masks using volume fractions: $g_{\text{agg}} = \sqrt{\text{AVF} / \text{FVF}}$ Stikov et al. (2015), then averaged across images within each region.

Columns with a red background are reference values from previous publications Stikov et al. (2015)Stikov et al. (2015).

from IPython.display import HTML

summary_rows = []
for sample in sorted(all_data.keys()):
    df = all_data[sample]
    ag = agg_gratio[sample]
    ref = PAPER_REF[sample]
    summary_rows.append({
        "Region": REGION_LABELS[sample],
        "n": len(df),
        "Fiber diam (\u03bcm)": f"{df['fiber_diam'].mean():.2f} \u00b1 {df['fiber_diam'].std():.2f}",
        "G-ratio (per-axon)": f"{df['gratio'].mean():.2f} \u00b1 {df['gratio'].std():.2f}",
        "G-ratio (aggregate)": f"{np.mean(ag):.2f} \u00b1 {np.std(ag):.2f}",
        "MRI g [S1]": ref["mri_g"],
        "Agg g (hist) [S1]": ref["agg_g_p1"],
        "Agg g [S2]": ref["agg_g_p2"],
        "Mean g [S2]": ref["mean_g_p2"],
    })

summary_df = pd.DataFrame(summary_rows)

paper_cols = ["MRI g [S1]", "Agg g (hist) [S1]", "Agg g [S2]", "Mean g [S2]"]
styled = summary_df.style.map(
    lambda _: "background-color: #fde0e0",
    subset=paper_cols
).hide(axis="index")

display(HTML(styled.to_html()))

Axon diameter vs Fiber diameter¶

A linear regression is fit to assess the linearity of the axon-fiber relationship. Use the slider to switch between CC regions.

samples_sorted = sorted(all_data.keys())
all_df = pd.concat(all_data.values())
x_max = all_df["axon_diam"].max() * 1.05
y_max = all_df["fiber_diam"].max() * 1.05

fig = go.Figure()

for i, sample in enumerate(samples_sorted):
    df = all_data[sample]
    x, y = df["axon_diam"].values, df["fiber_diam"].values
    visible = (i == 0)

    fig.add_trace(go.Scattergl(
        x=x, y=y, mode="markers",
        marker=dict(size=4, color="steelblue", opacity=0.3),
        showlegend=False, visible=visible,
        hovertemplate="Axon: %{x:.3f} \u03bcm<br>Fiber: %{y:.3f} \u03bcm<extra></extra>",
    ))

n_traces = len(samples_sorted)
steps = []
for i, sample in enumerate(samples_sorted):
    df = all_data[sample]
    vis = [False] * n_traces
    vis[i] = True
    step = dict(
        method="update",
        args=[
            {"visible": vis},
            {"title": f"{REGION_LABELS[sample]} (n = {len(df)})"},
        ],
        label=REGION_LABELS[sample],
    )
    steps.append(step)

fig.update_layout(
    sliders=[dict(active=0, steps=steps, currentvalue=dict(prefix="Region: "), pad=dict(t=50))],
    title=f"{REGION_LABELS[samples_sorted[0]]} (n = {len(all_data[samples_sorted[0]])})",
    xaxis_title="Axon diameter (\u03bcm)",
    yaxis_title="Fiber diameter (\u03bcm)",
    xaxis_range=[0, x_max], yaxis_range=[0, y_max],
    template=PLOT_TEMPLATE, height=600, width=800,
)
fig.show()

G-ratio vs Fiber diameter¶

The g-ratio is defined as the ratio of axon diameter to fiber diameter. A linear fit is applied to each region. Use the slider to switch between CC regions.

x_max_g = all_df["fiber_diam"].max() * 1.05
y_min_g = max(all_df["gratio"].min() * 0.95, 0)
y_max_g = min(all_df["gratio"].max() * 1.05, 1.0)

fig = go.Figure()

for i, sample in enumerate(samples_sorted):
    df = all_data[sample]
    x, y = df["fiber_diam"].values, df["gratio"].values
    visible = (i == 0)

    fig.add_trace(go.Scattergl(
        x=x, y=y, mode="markers",
        marker=dict(size=4, color="darkorange", opacity=0.3),
        showlegend=False, visible=visible,
        hovertemplate="Fiber: %{x:.3f} \u03bcm<br>G-ratio: %{y:.3f}<extra></extra>",
    ))

n_traces = len(samples_sorted)
steps = []
for i, sample in enumerate(samples_sorted):
    df = all_data[sample]
    vis = [False] * n_traces
    vis[i] = True
    step = dict(
        method="update",
        args=[
            {"visible": vis},
            {"title": f"{REGION_LABELS[sample]} (n = {len(df)})"},
        ],
        label=REGION_LABELS[sample],
    )
    steps.append(step)

fig.update_layout(
    sliders=[dict(active=0, steps=steps, currentvalue=dict(prefix="Region: "), pad=dict(t=50))],
    title=f"{REGION_LABELS[samples_sorted[0]]} (n = {len(all_data[samples_sorted[0]])})",
    xaxis_title="Fiber diameter (\u03bcm)",
    yaxis_title="G-ratio",
    xaxis_range=[0, x_max_g], yaxis_range=[y_min_g, y_max_g],
    template=PLOT_TEMPLATE, height=600, width=800,
)
fig.show()

References¶

Stikov, N., Campbell, J. S. W., Stroh, T., Lavelée, M., Frey, S., Novek, J., Nuara, S., Ho, M.-K., Bedell, B. J., Dougherty, R. F., Leppert, I. R., Boudreau, M., Narayanan, S., Duval, T., Cohen-Adad, J., Picard, P.-A., Gasecka, A., Côté, D., & Pike, G. B. (2015). In vivo histology of the myelin g-ratio with magnetic resonance imaging. NeuroImage, 118, 397–405. 10.1016/j.neuroimage.2015.05.023
Stikov, N., Campbell, J. S. W., Stroh, T., Lavelée, M., Frey, S., Novek, J., Nuara, S., Ho, M.-K., Bedell, B. J., Dougherty, R. F., Leppert, I. R., Boudreau, M., Narayanan, S., Duval, T., Cohen-Adad, J., Picard, P.-A., Gasecka, A., Côté, D., & Pike, G. B. (2015). Quantitative analysis of the myelin g-ratio from electron microscopy images of the macaque corpus callosum. Data in Brief, 4, 368–373. 10.1016/j.dib.2015.05.019