Parkinson’s disease in the spinal cord: an exploratory study to establish T2*w, MTR and diffusion-weighted imaging metric values

import plotly.express as px
from plotly.offline import plot
from IPython.core.display import HTML
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import pandas as pd
import statsmodels.formula.api as smf
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import base64
import plotly.io as pio
pio.renderers.default = "plotly_mimetype"

# Initialize the figure with subplots
fig = make_subplots(
    rows=1, cols=2, 
    vertical_spacing=0,  # Adjust vertical spacing (between rows)
    horizontal_spacing=0,  # Adjust horizontal spacing (between columns)
)

# Color palettes
palette_boxplots = ['steelblue', '#F0B0B0', 'lightcoral', '#B4464F']
palette_points = ['#00517F', '#B4464F', '#B4464F', '#5E000E']

# Load the data
df_CSA = pd.read_csv('../../data/parkinsons-spinalcord-mri-metrics/data/CSA.csv')
df_CSA = df_CSA[df_CSA['SpinalLevel'] == '2:05'] # Sort data for C2-C5 spinal levels combined 

# Define the groups
groups = ['CTRL', 'early', 'mid', 'adv']

# OLS analysis   
ols_model = smf.ols(formula='CSA ~ C(UPDRS_class_bis) + Age', data=df_CSA) # The "C()" here means that UPDRS_class_bis is a categorical variable
ols_results = ols_model.fit()
#print(f'OLS results for WM MTR in C2-C5: {ols_results.summary()}')

# Perform ANOVA test
anova_results = anova_lm(ols_results, typ=2)  # Type II ANOVA
#print(f'ANOVA results : {anova_results}')

# Adjusted R2
adjR2 = ols_results.rsquared_adj

# UPDRSIII_class_bis p-value
pvalue_UPDRS_class_bis = anova_results.loc['C(UPDRS_class_bis)', 'PR(>F)']

# Age p-value
pvalue_age = anova_results.loc['Age', 'PR(>F)']
formatted_pvalue_age = f"{pvalue_age:.4f}"

for i, group in enumerate(groups):
    box_color = palette_boxplots[i % len(palette_boxplots)]  
    points_color = palette_points[i % len(palette_points)] 
    
    fig.add_trace(go.Box(
        # Filter data for each group
        x=df_CSA['UPDRS_class_bis'][df_CSA['UPDRS_class_bis'] == group], 
        y=df_CSA['CSA'][df_CSA['UPDRS_class_bis'] == group],  

        # Show all points
        boxpoints='all',
        jitter=0.5,  # Jitter the points for better visibility
        whiskerwidth=0.8,  
        fillcolor=box_color,  # Color of the box
        marker_size=6,  # Marker size for points
        marker_color=points_color, # Color of the points
        line_width=2.5,  # Border width of the box
        line_color="black" ,  # Border color of the box
        pointpos=0, # Center the points with the box
        showlegend=False
    ), row=1, col=2
    )

# Determine if the p-value is significant (whether to add an asterisk on the plot or not)
red_asterisk = '<span style="color:red; font-size:30">*</span>' if pvalue_UPDRS_class_bis < 0.05 else ""
black_asterisk = '<span style="color:black; font-size:30">*</span>' if pvalue_age < 0.05 else ""

# Add annotations for p-values
fig.add_annotation(
x=0.98,  # Position of the annotation in x
y=0.98,  # Position of the annotation in y 
text=f"p-Group: {pvalue_UPDRS_class_bis:.4f}{red_asterisk}<br>p-Age: {pvalue_age:.4f}{black_asterisk}", 
showarrow=False,
font=dict(size=20, family="Arial", color="black"),
bgcolor="white",  # Background color of the annotation
align="right",
row=1,
col=2,
xref="x domain", 
yref="y domain"  
)

# Load static background image (from the "templates_for_figures" folder) and encode as base64
with open("../templates_for_figures/figure5_template.png", "rb") as image_file:
    encoded_image = base64.b64encode(image_file.read()).decode()

fig.update_layout(
    width=800, # Width of the overall figure
    height=450, # Height of the overall figure
    yaxis2_title="CSA (mm²)",
    yaxis2_title_font=dict(size=45, family="Arial", color="black", weight='bold'),
    images=[dict(
        source="data:image/png;base64," + encoded_image, # Add static background image
        x=-0.03,
        y=-0.16,
        xanchor="left",
        yanchor="bottom",
        sizex=1.3,
        sizey=1.3, 
        layer="below",  
    )],
)

fig.update_xaxes(
        ticktext=['HC', 'Early', 'Mid', 'Adv'], 
        tickvals=['CTRL', 'early', 'mid', 'adv'], 
        tickfont=dict(size=25, weight='bold', family="Arial"),
        row=1, col=2
    )

fig.update_yaxes(range=[60,100], tickfont=dict(size=27), row=1, col=2)

fig.show()