Show single subject data together with normative values

statistics/generate_figures_with_single_subject.py

@@ -0,0 +1,327 @@
+#
+# Functions to plot single subject morphometric metrics together with normative values computed from normative
+# database (spine-generic dataset in PAM50 space) per slice and vertebral levels
+#
+# Example usage:
+#       python generate_figures.py
+#       -path-HC $SCT_DIR/data/PAM50_normalized_metrics
+#       -participant-file $SCT_DIR/data/PAM50_normalized_metrics/participants.tsv
+#       -single-subject-sex M
+#
+# Author: Jan Valosek
+#
+
+import os
+import sys
+import re
+import argparse
+import numpy as np
+import pandas as pd
+import seaborn as sns
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+
+METRICS = ['MEAN(area)', 'MEAN(diameter_AP)', 'MEAN(diameter_RL)', 'MEAN(compression_ratio)', 'MEAN(eccentricity)',
+           'MEAN(solidity)']
+
+METRICS_DTYPE = {
+    'MEAN(diameter_AP)': 'float64',
+    'MEAN(area)': 'float64',
+    'MEAN(diameter_RL)': 'float64',
+    'MEAN(eccentricity)': 'float64',
+    'MEAN(solidity)': 'float64'
+}
+
+METRIC_TO_AXIS = {
+    'MEAN(diameter_AP)': 'AP Diameter [mm]',
+    'MEAN(area)': 'Cross-Sectional Area [mm²]',
+    'MEAN(diameter_RL)': 'Transverse Diameter [mm]',
+    'MEAN(eccentricity)': 'Eccentricity [a.u.]',
+    'MEAN(solidity)': 'Solidity [%]',
+    'MEAN(compression_ratio)': 'Compression Ratio [a.u.]',
+}
+
+LABELS_FONT_SIZE = 14
+TICKS_FONT_SIZE = 12
+
+AGE_DECADES = ['10-20', '21-30', '31-40', '41-50', '51-60']
+
+COLORS_SEX = {
+    'M': 'blue',
+    'F': 'red'
+    }
+
+SEX_TO_LEGEND = {
+    'M': 'males',
+    'F': 'females'
+    }
+
+
+def get_parser():
+    parser = argparse.ArgumentParser(
+        description="Plot single subject morphometric metrics together with normative values computed from normative "
+                    "database (spine-generic dataset in PAM50 space) per slice and vertebral levels ")
+    parser.add_argument('-path-HC', required=True, type=str,
+                        help="Path to data of normative dataset computed per slice. "
+                             "Example: '${SCT_DIR}/data/PAM50_normalized_metrics'")
+    parser.add_argument('-csv-single-subject', required=True, type=str,
+                        help="Path to the CSV file with single-subject morphometric metrics in the PAM50 space. "
+                             "The file can be generated using the 'sct_process_segmentation -perslice 1 "
+                             "-normalize-PAM50 1' command. "
+                             "Example: 'sub-001_T2w_metrics_PAM50.csv'")
+    parser.add_argument('-single-subject-sex', required=False, type=str, choices=['M', 'F'], default=None,
+                        help="Sex of the single subject. Options: 'M', 'F'. "
+                             "If used, the '-participant-file' flag must be also specified.")
+    parser.add_argument('-participant-file', required=False, type=str,
+                        help="Path to the spine-generic participants.tsv file (used to filter per sex). "
+                             "Example: '${SCT_DIR}/data/PAM50_normalized_metrics/participants.tsv'")
+    parser.add_argument('-path-out', required=False, type=str, default='figures',
+                        help="Output directory name. Default: figures.")
+
+    return parser
+
+
+def get_vert_indices(df):
+    """
+    Get indices of slices corresponding to mid-vertebrae
+    Args:
+        df (pd.dataFrame): dataframe with CSA values
+    Returns:
+        vert (pd.Series): vertebrae levels across slices
+        ind_vert (np.array): indices of slices corresponding to the beginning of each level (=intervertebral disc)
+        ind_vert_mid (np.array): indices of slices corresponding to mid-levels
+    """
+    # Get vert levels for one certain subject
+    vert = df[df['participant_id'] == 'sub-amu01']['VertLevel']
+    # Get indexes of where array changes value
+    ind_vert = vert.diff()[vert.diff() != 0].index.values
+    # Get the beginning of C1
+    ind_vert = np.append(ind_vert, vert.index.values[-1])
+    ind_vert_mid = []
+    # Get indexes of mid-vertebrae
+    for i in range(len(ind_vert)-1):
+        ind_vert_mid.append(int(ind_vert[i:i+2].mean()))
+
+    return vert, ind_vert, ind_vert_mid
+
+
+def fetch_subject_and_session(filename_path):
+    """
+    Get subject ID, session ID and filename from the input BIDS-compatible filename or file path
+    The function works both on absolute file path as well as filename
+    :param filename_path: input nifti filename (e.g., sub-001_ses-01_T1w.nii.gz) or file path
+    (e.g., /home/user/MRI/bids/derivatives/labels/sub-001/ses-01/anat/sub-001_ses-01_T1w.nii.gz
+    :return: subjectID: subject ID (e.g., sub-001)
+    """
+
+    _, filename = os.path.split(filename_path)              # Get just the filename (i.e., remove the path)
+    subject = re.search('sub-(.*?)[_/]', filename_path)
+    subjectID = subject.group(0)[:-1] if subject else ""    # [:-1] removes the last underscore or slash
+
+    session = re.search('ses-(.*?)[_/]', filename_path)     # [_/] means either underscore or slash
+    sessionID = session.group(0)[:-1] if session else ""    # [:-1] removes the last underscore or slash
+    # REGEX explanation
+    # \d - digit
+    # \d? - no or one occurrence of digit
+    # *? - match the previous element as few times as possible (zero or more times)
+
+    return subjectID, sessionID
+
+
+def create_lineplot(df, df_single_subject, sub_ses, number_of_subjects, path_out, sex=None):
+    """
+    Create lineplot for individual metrics per vertebral levels.
+    Note: we are plotting slices not levels to avoid averaging across levels.
+    Args:
+        df (pd.dataFrame): dataframe with metric values
+        df_single_subject (pd.dataFrame): dataframe with metric values for a single subject
+        sub_ses (str): subject and session ID
+        number_of_subjects (int): number of subjects in spine-generic dataset
+        path_out (str): path to output directory
+        sex (str): sex to filter spine-generic subjects on; possible options: 'M', 'F'
+    """
+
+    mpl.rcParams['font.family'] = 'Arial'
+
+    fig, axes = plt.subplots(2, 3, figsize=(20, 10))
+    axs = axes.ravel()
+
+    # Loop across metrics
+    for index, metric in enumerate(METRICS):
+        # Note: we are plotting slices not levels to avoid averaging across levels
+        if sex:
+            # Plot spine-generic multi-subject data for a given sex
+            sns.lineplot(ax=axs[index], x="Slice (I->S)", y=metric, data=df[df['sex'] == sex], errorbar='sd',
+                         linewidth=2, color=COLORS_SEX[sex],
+                         label=f'spine-generic {SEX_TO_LEGEND[sex]} (N = {number_of_subjects})')
+        else:
+            # Plot spine-generic multi-subject data
+            sns.lineplot(ax=axs[index], x="Slice (I->S)", y=metric, data=df, errorbar='sd', linewidth=2,
+                         label=f'spine-generic all subjects (N = {number_of_subjects})')
+        # Plot single subject data
+        sns.lineplot(ax=axs[index], x="Slice (I->S)", y=metric, data=df_single_subject, linewidth=2, color='black',
+                     label=sub_ses)
+
+        ymin, ymax = axs[index].get_ylim()
+
+        # Add legend
+        if index == 1:
+            axs[index].legend(loc='upper right', fontsize=TICKS_FONT_SIZE)
+        else:
+            axs[index].get_legend().remove()
+
+        # Add master title
+        plt.suptitle(f'Morphometric measures for {sub_ses} in the PAM50 template space', fontweight='bold',
+                     fontsize=LABELS_FONT_SIZE, y=0.92)
+
+        # Add labels
+        axs[index].set_ylabel(METRIC_TO_AXIS[metric], fontsize=LABELS_FONT_SIZE)
+        axs[index].set_xlabel('PAM50 Axial Slice #', fontsize=LABELS_FONT_SIZE)
+        # Increase xticks and yticks font size
+        axs[index].tick_params(axis='both', which='major', labelsize=TICKS_FONT_SIZE)
+
+        # Remove spines
+        axs[index].spines['right'].set_visible(False)
+        axs[index].spines['left'].set_visible(False)
+        axs[index].spines['top'].set_visible(False)
+        axs[index].spines['bottom'].set_visible(True)
+
+        # Get indices of slices corresponding vertebral levels
+        vert, ind_vert, ind_vert_mid = get_vert_indices(df)
+        # Insert a vertical line for each intervertebral disc
+        for idx, x in enumerate(ind_vert[1:-1]):
+            axs[index].axvline(df.loc[x, 'Slice (I->S)'], color='black', linestyle='--', alpha=0.5, zorder=0)
+
+        # Insert a text label for each vertebral level
+        for idx, x in enumerate(ind_vert_mid, 0):
+            # Deal with T1 label (C8 -> T1)
+            if vert[x] > 7:
+                level = 'T' + str(vert[x] - 7)
+                axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], ymin, level, horizontalalignment='center',
+                                verticalalignment='bottom', color='black', fontsize=TICKS_FONT_SIZE)
+            else:
+                level = 'C' + str(vert[x])
+                axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], ymin, level, horizontalalignment='center',
+                                verticalalignment='bottom', color='black', fontsize=TICKS_FONT_SIZE)
+
+        # Invert x-axis
+        axs[index].invert_xaxis()
+        # Add only horizontal grid lines
+        axs[index].yaxis.grid(True)
+        # Move grid to background (i.e. behind other elements)
+        axs[index].set_axisbelow(True)
+
+    # Save figure
+    filename = sub_ses + '_T2w_lineplot_PAM50.png'
+    path_filename = os.path.join(path_out, filename)
+    plt.savefig(path_filename, dpi=300, bbox_inches='tight')
+    print('Figure saved: ' + path_filename)
+
+
+def compute_cv(df, metric):
+    """
+    Compute coefficient of variation (CV) of a given metric.
+    Args:
+        df (pd.dataFrame): dataframe with CSA values
+        metric (str): column name of the dataframe to compute CV
+    Returns:
+        cv (float): coefficient of variation
+    """
+    cv = df[metric].std() / df[metric].mean()
+    cv = cv * 100
+    return cv
+
+
+def main():
+    parser = get_parser()
+    args = parser.parse_args()
+    path_HC = args.path_HC
+    path_single_subject = args.csv_single_subject
+    single_subject_sex = args.single_subject_sex
+    path_out_figures = os.path.abspath(args.path_out)
+
+    # Check that the arg -single-subject-sex is not used without -participant-file
+    if args.single_subject_sex and not args.participant_file:
+        print("ERROR: '-participant-file' flag is required when using the '-single-subject-sex' flag. Exiting...")
+        sys.exit(1)
+
+    # If the output folder directory is not present, then create it.
+    if not os.path.exists(path_out_figures):
+        os.makedirs(path_out_figures)
+    # Initialize pandas dataframe where data across all subjects will be stored
+    df = pd.DataFrame()
+    # Loop through .csv files of healthy controls
+    for file in os.listdir(path_HC):
+        if 'PAM50.csv' in file:
+            # Read csv file as pandas dataframe for given subject
+            df_subject = pd.read_csv(os.path.join(path_HC, file), dtype=METRICS_DTYPE)
+
+            # Compute compression ratio (CR) as MEAN(diameter_AP) / MEAN(diameter_RL)
+            df_subject['MEAN(compression_ratio)'] = df_subject['MEAN(diameter_AP)'] / df_subject['MEAN(diameter_RL)']
+
+            # Concatenate DataFrame objects
+            df = pd.concat([df, df_subject], axis=0, ignore_index=True)
+    # Get sub-id (e.g., sub-amu01) from Filename column and insert it as a new column called participant_id
+    # Subject ID is the first characters of the filename till slash
+    df.insert(0, 'participant_id', df['Filename'].str.split('/').str[0])
+
+    # If a participants.tsv file is provided, insert columns sex, age and manufacturer from df_participants into df
+    if args.participant_file:
+        df_participants = pd.read_csv(args.participant_file, sep='\t')
+        df = df.merge(df_participants[["age", "sex", "height", "weight", "manufacturer", "participant_id"]],
+                      on='participant_id')
+        # Recode age into age bins by 10 years (decades)
+        df['age'] = pd.cut(df['age'], bins=[10, 20, 30, 40, 50, 60], labels=AGE_DECADES)
+
+    df = df.dropna(axis=1, how='all')
+    df = df.dropna(axis=0, how='any').reset_index(drop=True)
+    # Keep only VertLevel from C1 to Th1
+    df = df[df['VertLevel'] <= 8]
+
+    df_spine_generic_min, df_spine_generic_max = df['Slice (I->S)'].min(), df['Slice (I->S)'].max()
+
+    # Multiply solidity by 100 to get percentage (sct_process_segmentation computes solidity in the interval 0-1)
+    df['MEAN(solidity)'] = df['MEAN(solidity)'] * 100
+
+    # Uncomment to save aggregated dataframe with metrics across all subjects as .csv file
+    #df.to_csv(os.path.join(path_out_csv, 'HC_metrics.csv'), index=False)
+
+    # Load single subject data
+    df_single_subject = pd.read_csv(path_single_subject, dtype=METRICS_DTYPE)
+    # Compute compression ratio (CR) as MEAN(diameter_AP) / MEAN(diameter_RL)
+    df_single_subject['MEAN(compression_ratio)'] = df_single_subject['MEAN(diameter_AP)'] / \
+                                                   df_single_subject['MEAN(diameter_RL)']
+    # Multiply solidity by 100 to get percentage (sct_process_segmentation computes solidity in the interval 0-1)
+    df_single_subject['MEAN(solidity)'] = df_single_subject['MEAN(solidity)'] * 100
+
+    # Keep only slices from C1 to Th1 to match the slices of the spine-generic normative values
+    df_single_subject = df_single_subject[(df_single_subject['Slice (I->S)'] >= df_spine_generic_min) &
+                                          (df_single_subject['Slice (I->S)'] <= df_spine_generic_max)]
+
+    # Check if df_single_subject is not empty, if so, print warning and exit
+    if df_single_subject.empty:
+        print('WARNING: No slices found in the range C1-Th1 in the single subject data. Exiting...')
+        sys.exit(1)
+
+    # Get subject ID from filename
+    subjectID, sessionID = fetch_subject_and_session(path_single_subject)
+
+    if sessionID:
+        sub_ses = subjectID + '_' + sessionID
+    else:
+        sub_ses = subjectID
+
+    if args.single_subject_sex:
+        number_of_subjects = str(len(df_participants[df_participants['sex'] == single_subject_sex]))
+        print(f"Number of {SEX_TO_LEGEND[single_subject_sex]}: {number_of_subjects}")
+    else:
+        number_of_subjects = str(len(df_participants))
+        print(f'Number of subjects: {number_of_subjects}')
+
+    # Create plots
+    create_lineplot(df, df_single_subject, sub_ses, number_of_subjects, path_out_figures, sex=single_subject_sex)
+
+
+if __name__ == '__main__':
+    main()