2024_Ruthe_SND/code/plotting_quantitative_dti_values.py

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

@author: arefk
"""

import os
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from statsmodels.stats.multicomp import MultiComparison
from itertools import combinations
import statsmodels
import scipy
from scipy.stats import dunnett

# Set up warning suppression for convergence
import warnings
warnings.simplefilter(action='ignore')

# Define user settings for the analysis
selected_dialation_amount = 3  # Set dilation amount
selected_qtype = "fa"  # Set Qtype value
selected_mask_names = ["L_cst_784_AMBA",
    "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",
    "CC_MOp-MOp_cut",
    "OT_och-lgn_lgncut"
]  # Set mask names to be used

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

# Define the PlotBoxplot function
def PlotBoxplot(filtered_df, qq, mm, dd, output_folder, sig):
    """
    Plot boxplots for filtered dataframe and save as images.

    Parameters:
    filtered_df (DataFrame): Filtered dataframe based on "Qtype" and "mask_name".
    qq (str): Qtype value.
    mm (str): Mask_name value.
    output_folder (str): Path to the output folder to save images.
    """
    # Set the color palette to grey
    sns.set_palette("rocket")
    
    # Set the figure size to 18 cm in width
    plt.figure(figsize=(18/2.54, 4))  # Convert 18 cm to inches
    
    # Create a boxplot
    sns.boxplot(data=filtered_df, hue="Group", y="Value", x="merged_timepoint",
                fliersize=3, flierprops={"marker": "o"}, palette="rocket")
    
    
    # Set title and labels
    mask_label = mask_name_mapping.get(mm, mm)
    plt.title(f'{qq}_{mask_label}_{dd}')
    plt.xlabel('Timepoint')
    plt.ylabel(qq + "(n.a)")
    
    # Show legend without title and frame
    plt.legend(title=None, frameon=False, bbox_to_anchor=(1.05, 1), loc='best')
    plt.tight_layout()
    
    # Save the plot as an image
    if sig:
        output_path = os.path.join(output_folder, f'{qq}_{mm}_{dd}_mixedlm_significant.png')
    else:
        output_path = os.path.join(fa_over_time_plots_dir, f'{mm}.png')
    plt.savefig(output_path, dpi=100)
    
    # Close the plot to avoid displaying it
    plt.show()
    plt.close()

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

# Get the parent directory of the code directory
parent_dir = os.path.dirname(code_dir)

# Define the path for the input CSV file
input_file_path = os.path.join(parent_dir, 'output', "Quantitative_outputs","onlyCST", 'Quantitative_results_from_dwi_processing_withOutWM_and_partialCST.csv')

# Define the path for the output folder to save plots
plots_output_dir = os.path.join(parent_dir, 'output', 'Figures', 'pythonFigs')
fa_over_time_plots_dir = os.path.join(parent_dir, 'output', 'Figures', 'fa_over_time_plots')
os.makedirs(fa_over_time_plots_dir, exist_ok=True)
os.makedirs(plots_output_dir, exist_ok=True)

# Load the CSV file into a dataframe
df = pd.read_csv(input_file_path)
#df = df[df["cluster ii"]=="High ii"]
df = df.drop_duplicates()
df = df[df["merged_timepoint"].isin([0,3,7,14,21, 28])]


# Create lists to store results and terminal output
results = []

# Iterate over the user-defined "dialation_amount", "Qtype", and "mask_name"
for dd in [selected_dialation_amount]:
    for qq in [selected_qtype]:
        for mm in selected_mask_names:
            # Filter the dataframe based on current "Qtype", "mask_name", and "dialation_amount"
            filtered_df = df[(df["Qtype"] == qq) & (df["mask_name"] == mm) & (df["dialation_amount"] == dd)]
            
            # Check if the filtered dataframe is not empty
            if not filtered_df.empty:
                # Plot the boxplot for the filtered dataframe
                # Assuming `sig` is a boolean that indicates statistical significance (set to False for now)
                sig = False
                PlotBoxplot(filtered_df, qq, mm, dd, plots_output_dir, sig)

# Print the list of results (currently empty as calculations are not included)
print(results)  # Replace this with actual results dataframe if available