added significance tests

code/plotting_quantitative_dti_values.py

@@ -32,24 +32,28 @@ selected_mask_names_with_dilations = {
     "CST_MOp-int-py_ipsilesional_CCcut": 2,
     "OT_och-lgn_lgncut": 2,
     "TC_DORsm-SSp_ll+SSp_ul_ipsilesional_END+higherstepsize_": 0,
-    "TC_DORsm-SSp+SSs_contralesional_END+CCcut": 0
+    "TC_DORsm-SSp+SSs_contralesional_END+CCcut": 0,
+    "L_cst_784_AMBA":0,
+    "L_cc_776_AMBA":0,
+    "L_VIS_669_AMBA":0
 }  # Set mask names and their specific dilation amounts for Qtype 'fa'
 
-# Define simple anatomical names for the masks
+# Define simple anatomical names for the masks (this is not correct)
 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)",
-    "CC_MOp-MOp_cut": "Corpus Callosum",
-    "OT_och-lgn_lgncut": "Optic"
+    "CRusdT_MOp-RN_ipsilesional_CCcut": "Rubropsisdnal (Ipsilesional)",
+    "CRuT_MOp-RN_contralesional_sdmirrored_CCcut": "Rubsdropsinal (Contralesional)",
+    "CReT_MOp-TRN_ipsilesionalsd_CCcut": "Reticulsdospinal (Ipsilesional)",
+    "CReT_MOp-TRN_contralesionsdal_CCcut": "Reticulosdspinal (Contralesional)",
+    "CST_MOp-int-py_ipsilesional_selfdrawnROA+CsdCcut": "Csdorticospinal (Ipsilesional)",
+    "CST_MOp-int-py_contralesional_selfdrawsdnROA+CCcut": "Csdorticospinal (Contralesional)",
+    "CC_MOp-MOp_cusdt": "Corpus Calsdlosum",
+    "OT_och-lgn_lsdgncut": "Optisdc"
 }
 
 # Define the PlotBoxplot function
 
-# Define the PlotBoxplot function
+import numpy as np
+from pandas.api.types import CategoricalDtype
 
 def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder):
     """
@@ -69,6 +73,14 @@ def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder)
     # Set font to Calibri
     plt.rcParams['font.family'] = 'Calibri'
 
+    # Define the categorical order for merged_timepoint
+    timepoint_order = ["0", "3", "7", "14", "21", "28"]
+    cat_type = CategoricalDtype(categories=timepoint_order, ordered=True)
+
+    # Convert merged_timepoint to categorical with the defined order to align with categorical boxplot positions
+    filtered_df_stroke["merged_timepoint"] = filtered_df_stroke["merged_timepoint"].astype(str).astype(cat_type)
+    filtered_df_sham["merged_timepoint"] = filtered_df_sham["merged_timepoint"].astype(str).astype(cat_type)
+
     # Create the figure with subplots
     fig, axes = plt.subplots(1, 2, figsize=(18 / 2.54, 8 / 2.54))  # Convert cm to inches
 
@@ -79,11 +91,7 @@ def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder)
         whiskerprops={'color': 'black', 'linewidth': 0.8},
         fliersize=0, linewidth=1.0, boxprops={'linewidth': 0.8}
     )
-    sns.stripplot(
-        data=filtered_df_stroke, hue="Group", y="Value", x="merged_timepoint",
-        ax=axes[0], dodge=True, palette=custom_colors,
-        alpha=1.0, size=2, jitter=True, marker='o', legend=False
-    )
+
     sns.despine(ax=axes[0], right=True, top=True)
     axes[0].set_title(f'Stroke - {qq}\n{mask_name_mapping.get(mm, mm)}\nd={dd}', fontsize=12)
     axes[0].set_xlabel('Timepoint', fontsize=12)
@@ -92,6 +100,17 @@ def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder)
     axes[0].tick_params(axis='y', labelsize=10)
     axes[0].legend_.remove()  # Remove the legend from the Stroke plot
 
+    # Plot median line for Stroke group using matplotlib to ensure alignment
+    median_stroke = filtered_df_stroke.groupby("merged_timepoint")["Value"].median().reset_index()
+    axes[0].plot(median_stroke["merged_timepoint"], median_stroke["Value"], color='black', linewidth=1.0, marker='o', markersize=3)
+
+    # Plot strip plot for Stroke group last to make sure dots are on top
+    sns.stripplot(
+        data=filtered_df_stroke, hue="Group", y="Value", x="merged_timepoint",
+        ax=axes[0], dodge=True, palette={"Stroke": '#8b0000'},  # Darker red color for Stroke dots
+        alpha=1.0, size=2, jitter=True, marker='o', legend=False
+    )
+
     # Plot Sham group
     sns.boxplot(
         data=filtered_df_sham, hue="Group", y="Value", x="merged_timepoint",
@@ -99,11 +118,7 @@ def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder)
         whiskerprops={'color': 'black', 'linewidth': 0.8},
         fliersize=0, linewidth=1.0, boxprops={'linewidth': 0.8}
     )
-    sns.stripplot(
-        data=filtered_df_sham, hue="Group", y="Value", x="merged_timepoint",
-        ax=axes[1], dodge=True, palette=custom_colors,
-        alpha=1.0, size=2, jitter=True, marker='o', legend=False
-    )
+
     sns.despine(ax=axes[1], right=True, top=True)
     axes[1].set_title(f'Sham - {qq}\n{mask_name_mapping.get(mm, mm)}\nd={dd}', fontsize=12)
     axes[1].set_xlabel('Timepoint', fontsize=12)
@@ -112,6 +127,17 @@ def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder)
     axes[1].tick_params(axis='y', labelsize=10)
     axes[1].legend_.remove()  # Remove the legend from the Sham plot
 
+    # Plot median line for Sham group using matplotlib to ensure alignment
+    median_sham = filtered_df_sham.groupby("merged_timepoint")["Value"].median().reset_index()
+    axes[1].plot(median_sham["merged_timepoint"], median_sham["Value"], color='black', linewidth=1.0, marker='o', markersize=3)
+
+    # Plot strip plot for Sham group last to make sure dots are on top
+    sns.stripplot(
+        data=filtered_df_sham, hue="Group", y="Value", x="merged_timepoint",
+        ax=axes[1], dodge=True, palette={"Sham": '#4f4f4f'},  # Darker gray color for Sham dots
+        alpha=1.0, size=2, jitter=True, marker='o', legend=False
+    )
+
     plt.tight_layout()
 
     # Save the plot as an image
@@ -125,6 +151,11 @@ def PlotBoxplot(filtered_df_stroke, filtered_df_sham, qq, mm, dd, output_folder)
     plt.show()
     plt.close()
 
+
+
+
+
+
 # Get the directory where the code file is located
 code_dir = os.path.dirname(os.path.abspath(__file__))
 

code/ttest_BL_P3_differences.py

@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Oct 28 12:21:08 2024
+
+@author: arefk
+"""
+
+import os
+import pandas as pd
+import numpy as np
+from scipy.stats import shapiro, ttest_ind, mannwhitneyu
+from tqdm import tqdm
+from concurrent.futures import ThreadPoolExecutor, as_completed
+
+# 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)
+
+# Load the CSV data
+input_file_path = os.path.join(parent_dir, 'output', "Quantitative_outputs",'Quantitative_results_from_dwi_processing_only_in_stroke_affected_slices.csv')
+df = pd.read_csv(input_file_path, low_memory=False)
+
+# Initialize an empty list to store results
+results = []
+
+# Get unique values of masks, qtypes, groups, and dilation amounts
+unique_masks = df['mask_name'].unique()
+unique_qtypes = df['Qtype'].unique()
+unique_groups = df['Group'].unique()
+unique_dilations = df['dialation_amount'].unique()
+
+# Prepare the combinations for parallel processing
+combinations = [(mask, qtype, group, dilation) for mask in unique_masks for qtype in unique_qtypes for group in unique_groups for dilation in unique_dilations]
+
+# Function to process each combination
+def process_combination(mask, qtype, group, dilation):
+    result = None
+    
+    # Filter data for timepoints 0 and 3 separately
+    df_timepoint_0 = df[(df['merged_timepoint'] == 0) &
+                        (df['Group'] == group) &
+                        (df['mask_name'] == mask) &
+                        (df['Qtype'] == qtype) &
+                        (df['dialation_amount'] == dilation)]
+
+    df_timepoint_3 = df[(df['merged_timepoint'] == 3) &
+                         (df['Group'] == group) &
+                         (df['mask_name'] == mask) &
+                         (df['Qtype'] == qtype) &
+                         (df['dialation_amount'] == dilation)]
+
+    # Drop NaN values for the 'Value' column
+    timepoint_0_values = df_timepoint_0['Value'].dropna()
+    timepoint_3_values = df_timepoint_3['Value'].dropna()
+
+    # Filter data after dropping NaN values to get subjects with non-null values
+    df_timepoint_0_filtered = df_timepoint_0[df_timepoint_0['Value'].notna()]
+    df_timepoint_3_filtered = df_timepoint_3[df_timepoint_3['Value'].notna()]
+
+    # Only proceed if there are more than 8 subjects in either timepoint after dropping NaNs
+    if len(df_timepoint_0_filtered['subjectID'].unique()) > 8 and len(df_timepoint_3_filtered['subjectID'].unique()) > 8:
+
+        # Check if we have enough values to perform statistical tests
+        if len(timepoint_0_values) > 0 and len(timepoint_3_values) > 0:
+            # Perform Shapiro-Wilk normality test
+            shapiro_timepoint_0_p = shapiro(timepoint_0_values)[1]
+            shapiro_timepoint_3_p = shapiro(timepoint_3_values)[1]
+
+            # Check if data is normally distributed
+            if shapiro_timepoint_0_p < 0.05 or shapiro_timepoint_3_p < 0.05:
+                # Use Mann-Whitney U test if data is not normally distributed
+                stat, p_value = mannwhitneyu(timepoint_0_values, timepoint_3_values, alternative='two-sided')
+            else:
+                # Use Welch's t-test if data is normally distributed
+                stat, p_value = ttest_ind(timepoint_0_values, timepoint_3_values, equal_var=False)
+
+            # Store the result
+            result = {
+                'mask_name': mask,
+                'Qtype': qtype,
+                'Group': group,
+                'dialation_amount': dilation,
+                'Pvalue': p_value
+            }
+
+    return result
+
+# Parallel processing using ThreadPoolExecutor with 6 workers
+with ThreadPoolExecutor(max_workers=6) as executor:
+    futures = {executor.submit(process_combination, mask, qtype, group, dilation): (mask, qtype, group, dilation) for mask, qtype, group, dilation in combinations}
+
+    # Iterate through completed tasks with progress bar
+    with tqdm(total=len(futures), desc="Processing combinations in parallel") as pbar:
+        for future in as_completed(futures):
+            combination = futures[future]
+            try:
+                result = future.result()
+                if result:
+                    results.append(result)
+            except Exception as e:
+                print(f"Error processing combination {combination}: {e}")
+            finally:
+                pbar.update(1)
+
+# Convert results to a DataFrame
+results_df = pd.DataFrame(results)
+
+# Define output path for the new CSV
+output_file_path = os.path.join(parent_dir, 'output', "Quantitative_outputs", 'Significance_timepoint_0_vs_3_only_in_stroke_slices.csv')
+
+# Save results to CSV
+results_df.to_csv(output_file_path, index=False)
+
+print(f"Significance analysis results saved to {output_file_path}")

code/ttest_BL_P7_differences.py

@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Oct 28 12:21:08 2024
+
+@author: arefk
+"""
+
+import os
+import pandas as pd
+import numpy as np
+from scipy.stats import shapiro, ttest_ind, mannwhitneyu
+from tqdm import tqdm
+from concurrent.futures import ThreadPoolExecutor, as_completed
+
+# 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)
+
+# Load the CSV data
+input_file_path = os.path.join(parent_dir, 'output', "Quantitative_outputs",'Quantitative_results_from_dwi_processing_only_in_stroke_affected_slices.csv')
+df = pd.read_csv(input_file_path, low_memory=False)
+
+# Initialize an empty list to store results
+results_p0_p7 = []
+
+# Get unique values of masks, qtypes, groups, and dilation amounts
+unique_masks = df['mask_name'].unique()
+unique_qtypes = df['Qtype'].unique()
+unique_groups = df['Group'].unique()
+unique_dilations = df['dialation_amount'].unique()
+
+# Prepare the combinations for parallel processing
+combinations = [(mask, qtype, group, dilation) for mask in unique_masks for qtype in unique_qtypes for group in unique_groups for dilation in unique_dilations]
+
+# Function to process each combination for timepoints 0 and 7
+def process_combination_p0_p7(mask, qtype, group, dilation):
+    result = None
+    
+    # Filter data for timepoints 0 and 7 separately
+    df_timepoint_0 = df[(df['merged_timepoint'] == 0) &
+                        (df['Group'] == group) &
+                        (df['mask_name'] == mask) &
+                        (df['Qtype'] == qtype) &
+                        (df['dialation_amount'] == dilation)]
+
+    df_timepoint_7 = df[(df['merged_timepoint'] == 7) &
+                         (df['Group'] == group) &
+                         (df['mask_name'] == mask) &
+                         (df['Qtype'] == qtype) &
+                         (df['dialation_amount'] == dilation)]
+
+    # Drop NaN values for the 'Value' column
+    timepoint_0_values = df_timepoint_0['Value'].dropna()
+    timepoint_7_values = df_timepoint_7['Value'].dropna()
+
+    # Filter data after dropping NaN values to get subjects with non-null values
+    df_timepoint_0_filtered = df_timepoint_0[df_timepoint_0['Value'].notna()]
+    df_timepoint_7_filtered = df_timepoint_7[df_timepoint_7['Value'].notna()]
+
+    # Only proceed if there are more than 8 subjects in either timepoint after dropping NaNs
+    if len(df_timepoint_0_filtered['subjectID'].unique()) > 8 and len(df_timepoint_7_filtered['subjectID'].unique()) > 8:
+
+        # Check if we have enough values to perform statistical tests
+        if len(timepoint_0_values) > 0 and len(timepoint_7_values) > 0:
+            # Perform Shapiro-Wilk normality test
+            shapiro_timepoint_0_p = shapiro(timepoint_0_values)[1]
+            shapiro_timepoint_7_p = shapiro(timepoint_7_values)[1]
+
+            # Check if data is normally distributed
+            if shapiro_timepoint_0_p < 0.05 or shapiro_timepoint_7_p < 0.05:
+                # Use Mann-Whitney U test if data is not normally distributed
+                stat, p_value = mannwhitneyu(timepoint_0_values, timepoint_7_values, alternative='two-sided')
+            else:
+                # Use Welch's t-test if data is normally distributed
+                stat, p_value = ttest_ind(timepoint_0_values, timepoint_7_values, equal_var=False)
+
+            # Store the result
+            result = {
+                'mask_name': mask,
+                'Qtype': qtype,
+                'Group': group,
+                'dialation_amount': dilation,
+                'Pvalue': p_value
+            }
+
+    return result
+
+# Parallel processing using ThreadPoolExecutor with 6 workers
+with ThreadPoolExecutor(max_workers=6) as executor:
+    futures_p0_p7 = {executor.submit(process_combination_p0_p7, mask, qtype, group, dilation): (mask, qtype, group, dilation) for mask, qtype, group, dilation in combinations}
+
+    # Iterate through completed tasks with progress bar
+    with tqdm(total=len(futures_p0_p7), desc="Processing combinations for timepoints 0 vs 7 in parallel") as pbar:
+        for future in as_completed(futures_p0_p7):
+            combination = futures_p0_p7[future]
+            try:
+                result = future.result()
+                if result:
+                    results_p0_p7.append(result)
+            except Exception as e:
+                print(f"Error processing combination {combination}: {e}")
+            finally:
+                pbar.update(1)
+
+# Convert results to a DataFrame
+results_p0_p7_df = pd.DataFrame(results_p0_p7)
+
+# Define output path for the new CSV
+output_file_path_p0_p7 = os.path.join(parent_dir, 'output', "Quantitative_outputs", 'Significance_timepoint_0_vs_7_only_in_stroke_slices.csv')
+
+# Save results to CSV
+results_p0_p7_df.to_csv(output_file_path_p0_p7, index=False)
+
+print(f"Significance analysis results saved to {output_file_path_p0_p7}")

code/ttest_P3_P7_differences.py

@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Oct 28 12:21:08 2024
+
+@author: arefk
+"""
+
+import os
+import pandas as pd
+import numpy as np
+from scipy.stats import shapiro, ttest_ind, mannwhitneyu
+from tqdm import tqdm
+from concurrent.futures import ThreadPoolExecutor, as_completed
+
+# 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)
+
+# Load the CSV data
+input_file_path = os.path.join(parent_dir, 'output', "Quantitative_outputs",'Quantitative_results_from_dwi_processing_only_in_stroke_affected_slices.csv')
+df = pd.read_csv(input_file_path, low_memory=False)
+
+# Initialize an empty list to store results
+results_p3_p7 = []
+
+# Get unique values of masks, qtypes, groups, and dilation amounts
+unique_masks = df['mask_name'].unique()
+unique_qtypes = df['Qtype'].unique()
+unique_groups = df['Group'].unique()
+unique_dilations = df['dialation_amount'].unique()
+
+# Prepare the combinations for parallel processing
+combinations = [(mask, qtype, group, dilation) for mask in unique_masks for qtype in unique_qtypes for group in unique_groups for dilation in unique_dilations]
+
+# Function to process each combination for timepoints 3 and 7
+def process_combination_p3_p7(mask, qtype, group, dilation):
+    result = None
+    
+    # Filter data for timepoints 3 and 7 separately
+    df_timepoint_3 = df[(df['merged_timepoint'] == 3) &
+                        (df['Group'] == group) &
+                        (df['mask_name'] == mask) &
+                        (df['Qtype'] == qtype) &
+                        (df['dialation_amount'] == dilation)]
+
+    df_timepoint_7 = df[(df['merged_timepoint'] == 7) &
+                         (df['Group'] == group) &
+                         (df['mask_name'] == mask) &
+                         (df['Qtype'] == qtype) &
+                         (df['dialation_amount'] == dilation)]
+
+    # Drop NaN values for the 'Value' column
+    timepoint_3_values = df_timepoint_3['Value'].dropna()
+    timepoint_7_values = df_timepoint_7['Value'].dropna()
+
+    # Filter data after dropping NaN values to get subjects with non-null values
+    df_timepoint_3_filtered = df_timepoint_3[df_timepoint_3['Value'].notna()]
+    df_timepoint_7_filtered = df_timepoint_7[df_timepoint_7['Value'].notna()]
+
+    # Only proceed if there are more than 8 subjects in either timepoint after dropping NaNs
+    if len(df_timepoint_3_filtered['subjectID'].unique()) > 8 and len(df_timepoint_7_filtered['subjectID'].unique()) > 8:
+
+        # Check if we have enough values to perform statistical tests
+        if len(timepoint_3_values) > 0 and len(timepoint_7_values) > 0:
+            # Perform Shapiro-Wilk normality test
+            shapiro_timepoint_3_p = shapiro(timepoint_3_values)[1]
+            shapiro_timepoint_7_p = shapiro(timepoint_7_values)[1]
+
+            # Check if data is normally distributed
+            if shapiro_timepoint_3_p < 0.05 or shapiro_timepoint_7_p < 0.05:
+                # Use Mann-Whitney U test if data is not normally distributed
+                stat, p_value = mannwhitneyu(timepoint_3_values, timepoint_7_values, alternative='two-sided')
+            else:
+                # Use Welch's t-test if data is normally distributed
+                stat, p_value = ttest_ind(timepoint_3_values, timepoint_7_values, equal_var=False)
+
+            # Store the result
+            result = {
+                'mask_name': mask,
+                'Qtype': qtype,
+                'Group': group,
+                'dialation_amount': dilation,
+                'Pvalue': p_value
+            }
+
+    return result
+
+# Parallel processing using ThreadPoolExecutor with 6 workers
+with ThreadPoolExecutor(max_workers=6) as executor:
+    futures_p3_p7 = {executor.submit(process_combination_p3_p7, mask, qtype, group, dilation): (mask, qtype, group, dilation) for mask, qtype, group, dilation in combinations}
+
+    # Iterate through completed tasks with progress bar
+    with tqdm(total=len(futures_p3_p7), desc="Processing combinations for timepoints 3 vs 7 in parallel") as pbar:
+        for future in as_completed(futures_p3_p7):
+            combination = futures_p3_p7[future]
+            try:
+                result = future.result()
+                if result:
+                    results_p3_p7.append(result)
+            except Exception as e:
+                print(f"Error processing combination {combination}: {e}")
+            finally:
+                pbar.update(1)
+
+# Convert results to a DataFrame
+results_p3_p7_df = pd.DataFrame(results_p3_p7)
+
+# Define output path for the new CSV
+output_file_path_p3_p7 = os.path.join(parent_dir, 'output', "Quantitative_outputs", 'Significance_timepoint_3_vs_7_only_in_stroke_slices.csv')
+
+# Save results to CSV
+results_p3_p7_df.to_csv(output_file_path_p3_p7, index=False)
+
+print(f"Significance analysis results saved to {output_file_path_p3_p7}")