2024_Ruthe_SND/code/vizualization_asymetry_index_fa.py

# -*- coding: utf-8 -*-
"""
Created on Tue Nov  5 14:25:19 2024

@author: arefks
"""

# -*- coding: utf-8 -*-
"""
Created on Mon Nov  4 17:54:46 2024

@author: arefk
"""

import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import pearsonr, spearmanr, shapiro
from statsmodels.stats.multitest import multipletests

# Get the directory where the code file is located
code_dir = os.path.dirname(os.path.abspath(__file__))
parent_dir = os.path.dirname(code_dir)

# Make sure these variables are set according to your selected settings
selected_groups = ["Stroke"]  # Example: Stroke group vs. Control group
selected_qtype = "fa"  # Example: fa Qtype
selected_mask_names = [
    "CRuT_MOp-RN_ipsilesional_CCcut",
    "CRuT_MOp-RN_contralesional_mirrored_CCcut",
    "CReT_MOp-TRN_ipsilesional_CCcut",
    "CReT_MOp-TRN_contralesional_CCcut",
    "CST_MOp-int-py_ipsilesional_selfdrawnROA+CCcut",
    "CST_MOp-int-py_contralesional_selfdrawnROA+CCcut",
    "TC_DORsm-SSp_ll+SSp_ul_ipsilesional_END+higherstepsize_",
    "TC_DORsm-SSp_ll+SSp_ul_contralesional_END+higherstepsize"
]  # Updated list of specific masks  # Updated list of specific masks
selected_timepoints = [3]  # Example: multiple selected timepoints
selected_dialation_amount = 3  # Example: dilation value for the time point

# Define simple anatomical names for the masks
mask_name_mapping = {
    "CC_MOp-MOp_cut": "Corpus Callosum",
    "CRuT_MOp-RN_ipsilesional_CCcut": "Rubropsinal (Ipsilesional)",
    "CRuT_MOp-RN_contralesional_mirrored_CCcut": "Rubropsinal (Contralesional)",
    "TC_DORsm-SSp_ll+SSp_ul_ipsilesional_END+higherstepsize_": "Thalamocortical (Ipsilesional)",
    "TC_DORsm-SSp_ll+SSp_ul_contralesional_END+higherstepsize": "Thalamocortical (Contralesional)",
    "CReT_MOp-TRN_contralesional_CCcut": "Reticulospinal (Contralesional)",
    "CReT_MOp-TRN_ipsilesional_CCcut": "Reticulospinal (Ipsilesional)",
    "CST_MOp-int-py_contralesional_selfdrawnROA+CCcut": "Corticospinal (Contralesional)",
    "CST_MOp-int-py_ipsilesional_selfdrawnROA+CCcut": "Corticospinal (Ipsilesional)",
    "OT_och-lgn_lgncut": "Optic"
}

# Load the CSV data
input_file_path = os.path.join(parent_dir, 'output', "Quantitative_outputs", 'Quantitative_results_from_dwi_processing_merged_with_behavior_data.csv')
df = pd.read_csv(input_file_path, low_memory=False)
# Remove duplicate rows
df = df.drop_duplicates()

# Define output paths
figures_folder = os.path.join(parent_dir, 'output', 'Figures')
pythonFigs_folder = os.path.join(figures_folder, 'pythonFigs')

# Create directories if they do not exist
os.makedirs(pythonFigs_folder, exist_ok=True)

# Store all p-values for multiple test correction
all_p_values = []

# Create the figure for subplots
fig, axes = plt.subplots(4, 1, figsize=(14 / 2.54, 20 / 2.54))  # 1 column, 4 rows  # 2 columns, 5 rows
axes = axes.flatten() if isinstance(axes, np.ndarray) else [axes]  # Flatten the axes array for easy iteration

# File to store all correlation results
correlation_text_path = os.path.join(pythonFigs_folder, 'fa_behavior_correlation_asymmetry_index.txt')
with open(correlation_text_path, 'w') as f:
    f.write('Asymmetry Index Correlation Results for All Masks\n')
    f.write('=' * 60 + '\n')

# Iterate over each pair of ipsilesional and contralesional mask names
for idx in range(0, len(selected_mask_names), 2):
    ipsilesional_mask_name = selected_mask_names[idx]
    contralesional_mask_name = selected_mask_names[idx + 1]

    for selected_timepoint in selected_timepoints:
        # Filter the DataFrame based on selected settings for each group
        df_filtered_ipsi_list = []
        df_filtered_contra_list = []
        for group in selected_groups:
            df_filtered_ipsi = df[(df['Group'] == group) &
                                  (df['Qtype'] == selected_qtype) &
                                  (df['mask_name'] == ipsilesional_mask_name) &
                                  (df['dialation_amount'] == selected_dialation_amount) &
                                  (df['merged_timepoint'] == selected_timepoint)]
            df_filtered_contra = df[(df['Group'] == group) &
                                    (df['Qtype'] == selected_qtype) &
                                    (df['mask_name'] == contralesional_mask_name) &
                                    (df['dialation_amount'] == selected_dialation_amount) &
                                    (df['merged_timepoint'] == selected_timepoint)]
            df_filtered_ipsi_list.append(df_filtered_ipsi)
            df_filtered_contra_list.append(df_filtered_contra)

        # Combine the filtered DataFrames
        df_filtered_ipsi = pd.concat(df_filtered_ipsi_list, ignore_index=True)
        df_filtered_contra = pd.concat(df_filtered_contra_list, ignore_index=True)

        # Remove NaNs and ensure we have matching data for asymmetry index calculation
        df_filtered_ipsi = df_filtered_ipsi.dropna(subset=['Value', 'DeficitScore'])
        df_filtered_contra = df_filtered_contra.dropna(subset=['Value', 'DeficitScore'])

        # Ensure that we have matching rows for both ipsilesional and contralesional masks
        if not df_filtered_ipsi.empty and not df_filtered_contra.empty:
            # Calculate the Asymmetry Index (AI)
            df_filtered_ipsi['AsymmetryIndex'] = (
                (df_filtered_contra['Value'].values - df_filtered_ipsi['Value'].values) /
                (0.5 * (df_filtered_contra['Value'].values + df_filtered_ipsi['Value'].values))
            )

            # Plot the data if available
            y_min = df_filtered_ipsi['DeficitScore'].min() - 0.1 * df_filtered_ipsi['DeficitScore'].max()
            y_max = df_filtered_ipsi['DeficitScore'].max() + 0.3 * df_filtered_ipsi['DeficitScore'].max()

            ax = axes[idx // 2]  # Select the appropriate subplot axis

            group_p_values = []
            correlation_results = []

            for group in selected_groups:
                group_data = df_filtered_ipsi[df_filtered_ipsi['Group'] == group]

                if not group_data.empty:
                    # Shapiro-Wilk test for normality
                    shapiro_statValue, shapiro_pvalueAI = shapiro(group_data['AsymmetryIndex'])
                    shapiro_statScore, shapiro_pvalueBehavior = shapiro(group_data['DeficitScore'])

                    # Calculate Pearson or Spearman based on normality test
                    if shapiro_pvalueAI < 0.05 or shapiro_pvalueBehavior < 0.05:
                        # Use Spearman if data is not normally distributed
                        correlation_coefficient, p_value = spearmanr(group_data['AsymmetryIndex'], group_data['DeficitScore'])
                    else:
                        # Use Pearson if data is normally distributed
                        correlation_coefficient, p_value = pearsonr(group_data['AsymmetryIndex'], group_data['DeficitScore'])

                    # Store p-value for multiple testing correction
                    group_p_values.append(p_value)
                    all_p_values.append(p_value)

                    # Store correlation result
                    correlation_results.append((group, correlation_coefficient, p_value))

            # Apply multiple test correction for this group
            corrected_p_values = multipletests(group_p_values, method='fdr_bh')[1]
            corrected_p_index = 0

            for group, correlation_coefficient, p_value in correlation_results:
                # Retrieve corrected p-value
                corrected_p_value = corrected_p_values[corrected_p_index]
                corrected_p_index += 1

                # Determine significance level
                pval_text = ''
                if corrected_p_value < 0.001:
                    pval_text = '***'
                elif corrected_p_value < 0.01:
                    pval_text = '**'
                elif corrected_p_value < 0.05:
                    pval_text = '*'
                else:
                    pval_text = ' '

                # Plot the regression
                sns.regplot(
                    x='AsymmetryIndex', y='DeficitScore', data=group_data, ci=95, ax=ax,
                    line_kws={'color': 'lightcoral', 'alpha': 0.6},  # Regression line in light red with higher transparency
                    scatter_kws={'color': 'red', 'alpha': 0.6, 's': 10},  # Scatter points in red
                    label=f'{group} (R={correlation_coefficient:.2f} {pval_text})'
                )
                # Add text annotation for R and significance level
                ax.text(0.05, 0.95, f'R={correlation_coefficient:.2f}{pval_text}', transform=ax.transAxes,
                        fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0))

                # Save correlation results to the text file
                with open(correlation_text_path, 'a') as f:
                    f.write(f'Mask Pair: {ipsilesional_mask_name} & {contralesional_mask_name} | Group: {group}\n')
                    f.write(f'Correlation Coefficient (uncorrected): {correlation_coefficient:.4f}\n')
                    f.write(f'P-Value (uncorrected): {p_value:.4e}\n')
                    f.write(f'P-Value (corrected): {corrected_p_value:.4e}\n')
                    f.write('-' * 40 + '\n')

            # Set consistent y-limits
            ax.set_ylim([y_min, y_max])

            # Set labels and customize font settings
            anatomical_name = mask_name_mapping.get(ipsilesional_mask_name, ipsilesional_mask_name)
            ax.set_ylabel('Deficit Score [a.u]', fontsize=12, fontname='Calibri')
            ax.set_xlabel(f'Asymmetry Index ({anatomical_name})', fontsize=12, fontname='Calibri')

            # Remove gridlines
            ax.grid(False)
            sns.despine(ax=ax)  # Remove the top and right spines for cleaner look

        else:
            print(f"No data available for the selected settings: Groups={selected_groups}, Qtype={selected_qtype}, Mask Pair=({ipsilesional_mask_name}, {contralesional_mask_name}), Dilation value={selected_dialation_amount}, Timepoint={selected_timepoint}.")

# Adjust layout
plt.tight_layout()

# Save the complete figure with all subplots
fig_filename = 'fa_behavior_correlation_asymmetry_index_all_masks_timepoint_56_groups_comparison.svg'
plt.savefig(os.path.join(pythonFigs_folder, fig_filename), dpi=300, bbox_inches='tight', format="svg", transparent=True)

# Show the figure
plt.show()