From 337a5c702551d67f1adcf344aaa0e862521ef25d Mon Sep 17 00:00:00 2001 From: Mathieu Boudreau Date: Wed, 7 Feb 2024 12:18:01 -0400 Subject: [PATCH] Remove mis-included file --- content/data_processing-2.ipynb | 989 -------------------------------- 1 file changed, 989 deletions(-) delete mode 100644 content/data_processing-2.ipynb diff --git a/content/data_processing-2.ipynb b/content/data_processing-2.ipynb deleted file mode 100644 index d176545..0000000 --- a/content/data_processing-2.ipynb +++ /dev/null @@ -1,989 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "459f60d5", - "metadata": {}, - "source": [ - "# Analysis code for the paper \"RF shimming in the cervical spinal cord at 7T\"\n", - "\n", - "## Data\n", - "\n", - "The data can be downloaded from: https://openneuro.org/datasets/ds004906\n", - "\n", - "The structure of the input dataset is as follows (JSON sidecars are not listed for clarity):\n", - "~~~\n", - "ds004906\n", - "├── CHANGES\n", - "├── README\n", - "├── dataset_description.json\n", - "├── participants.json\n", - "├── participants.tsv\n", - "├── sub-01\n", - "│   ├── anat\n", - "│   │   ├── sub-01_acq-CP_T1w.nii.gz\n", - "│   │   ├── sub-01_acq-CP_T2starw.nii.gz\n", - "│   │   ├── sub-01_acq-CoV_T1w.nii.gz\n", - "│   │   ├── sub-01_acq-CoV_T2starw.nii.gz\n", - "│   │   ├── sub-01_acq-SAReff_T2starw.nii.gz\n", - "│   │   ├── sub-01_acq-patient_T2starw.nii.gz\n", - "│   │   ├── sub-01_acq-phase_T2starw.nii.gz\n", - "│   │   ├── sub-01_acq-target_T2starw.nii.gz\n", - "│   │   ├── sub-01_acq-volume_T2starw.nii.gz\n", - "│   └── fmap\n", - "│   ├── sub-01_acq-anatCP_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-anatCoV_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-anatSAReff_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-anatpatient_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-anatphase_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-anattarget_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-anatvolume_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-fampCP_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-fampCoV_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-fampSAReff_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-famppatient_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-fampphase_TB1TFL.nii.gz\n", - "│   ├── sub-01_acq-famptarget_TB1TFL.nii.gz\n", - "│   └── sub-01_acq-fampvolume_TB1TFL.nii.gz\n", - "├── sub-02\n", - "├── sub-03\n", - "├── sub-04\n", - "└── sub-05\n", - "~~~\n", - "\n", - "\n", - "## Overview of processing pipeline\n", - "\n", - "For each subject:\n", - "\n", - "- Process anat/T2starw (GRE)\n", - " - Segment the spinal cord (SC)\n", - " - Label vertebral levels using existing manual disc labels\n", - " - Create a mask of the cerebrospinal fluid (CSF)\n", - " - Extract the SC/CSF magnitude signal to assess the stability of the flip angle across shim methods\n", - "- Process fmap/TFL (flip angle maps)\n", - " - Register each B1 map (CP, CoV, etc.) to the GRE scan\n", - " - Apply the computed warping field to bring the segmentation and vertebral levels to the B1 map\n", - " - Convert the B1 map to nT/V units\n", - " - Extract the B1 map value within the SC\n", - "\n", - ">Slow processes are indicated with the emoji ⏳" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4f78c5ea", - "metadata": {}, - "outputs": [], - "source": [ - "# Necessary imports\n", - "\n", - "from datetime import datetime, timedelta\n", - "import os\n", - "import re\n", - "import json\n", - "import subprocess\n", - "import glob\n", - "import itertools\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import nibabel as nib\n", - "import pandas as pd\n", - "from scipy.interpolate import interp1d\n", - "from scipy.ndimage import uniform_filter1d\n", - "from scipy.stats import f_oneway, ttest_rel\n", - "import shutil\n", - "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n", - "from statsmodels.stats.anova import AnovaRM\n", - "import statsmodels.api as sm" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bc5a2ec0", - "metadata": {}, - "outputs": [], - "source": [ - "# Start timer\n", - "start_time = datetime.now()\n", - "\n", - "# Check where we are\n", - "!pwd" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5efab0b8", - "metadata": {}, - "outputs": [], - "source": [ - "# ⚠️ No need to run this cell if you run this notebook locally and already have these dependencies installed.\n", - "\n", - "# Install packages on system\n", - "!sudo apt-get update\n", - "!sudo apt-get install git-annex\n", - "\n", - "# Install Python libaries\n", - "!wget -O requirements.txt https://raw.githubusercontent.com/shimming-toolbox/rf-shimming-7t/main/requirements.txt\n", - "!pip install -r requirements.txt\n", - "\n", - "# Install SCT ⏳\n", - "!git clone --depth 1 https://github.com/spinalcordtoolbox/spinalcordtoolbox.git\n", - "!yes | spinalcordtoolbox/install_sct\n", - "os.environ['PATH'] += f\":/content/spinalcordtoolbox/bin\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3393a694", - "metadata": {}, - "outputs": [], - "source": [ - "# Download data and define path variables\n", - "\n", - "!datalad install https://github.com/OpenNeuroDatasets/ds004906.git\n", - "os.chdir(\"ds004906\")\n", - "!datalad get . # uncomment for production\n", - "# !datalad get sub-01/ # debugging\n", - "# Get derivatives containing manual labels\n", - "!datalad get derivatives" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6fe120e8", - "metadata": {}, - "outputs": [], - "source": [ - "# Define useful variables\n", - "\n", - "path_data = os.getcwd()\n", - "print(f\"path_data: {path_data}\")\n", - "path_labels = os.path.join(path_data, \"derivatives\", \"labels\")\n", - "path_qc = os.path.join(path_data, \"qc\")\n", - "shim_modes = [\"CP\", \"patient\", \"volume\", \"phase\", \"CoV\", \"target\", \"SAReff\"]\n", - "shim_modes_MPRAGE = [\"CP\", \"CoV\"] # TODO: change variable name PEP8\n", - "# shim_modes = [\"CP\", \"CoV\"] # debugging\n", - "print(f\"shim_modes: {shim_modes}\")\n", - "subjects = sorted(glob.glob(\"sub-*\"))\n", - "print(f\"subjects: {subjects}\")\n", - "\n", - "# Create output folder\n", - "path_results = os.path.join(path_data, 'derivatives', 'results')\n", - "os.makedirs(path_results, exist_ok=True)" - ] - }, - { - "cell_type": "markdown", - "id": "b0c4e225", - "metadata": {}, - "source": [ - "## Process anat/T2starw (GRE)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e95d87a4-3194-4369-b6cb-933ad1d5cfa3", - "metadata": {}, - "outputs": [], - "source": [ - "# Run segmentation on GRE scan\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " # The \"CoV\" RF shimming scenario was chosen as the segmentation baseline due to the more homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance in the C7-T2 region\n", - " fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-SC_seg.nii.gz\")\n", - " if os.path.exists(fname_manual_seg):\n", - " # Manual segmentation already exists. Copy it to local folder\n", - " print(f\"{subject}: Manual segmentation found\\n\")\n", - " shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_seg.nii.gz\")\n", - " # Generate QC report to make sure the manual segmentation is correct\n", - " !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n", - " else:\n", - " # Manual segmentation does not exist. Run automatic segmentation.\n", - " print(f\"{subject}: Manual segmentation not found\")\n", - " !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2s -qc {path_qc}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "eb9a9267", - "metadata": {}, - "outputs": [], - "source": [ - "# Copy CSF masks from the derivatives folder\n", - "\n", - "# For more details about how these masks were created, see: https://github.com/shimming-toolbox/rf-shimming-7t/issues/67\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n", - " if os.path.exists(fname_manual_seg):\n", - " # Manual segmentation already exists. Copy it to local folder\n", - " print(f\"{subject}: Manual segmentation found\\n\")\n", - " shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n", - " # Generate QC report to make sure the manual segmentation is correct\n", - " !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n", - " else:\n", - " # Manual segmentation does not exist. Panic!\n", - " print(f\"{subject}: Manual segmentation not found 😱\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "862dc562", - "metadata": {}, - "outputs": [], - "source": [ - "# Crop GRE scan for faster processing and better registration results\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " !sct_crop_image -i {subject}_acq-CoV_T2starw.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop.nii.gz\n", - " !sct_crop_image -i {subject}_acq-CoV_T2starw_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_seg.nii.gz\n", - " !sct_crop_image -i {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "87030540", - "metadata": {}, - "outputs": [], - "source": [ - "# Label vertebrae on GRE scan\n", - "\n", - "# Given the low resolution of the GRE scan, the automatic detection of C2-C3 disc is unreliable. Therefore we need to use the manual disc labels that are part of the dataset.\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " fname_label_discs = os.path.join(path_data, \"derivatives\", \"labels\", subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-discs_dseg.nii.gz\")\n", - " !sct_label_utils -i {subject}_acq-CoV_T2starw_crop_seg.nii.gz -disc {fname_label_discs} -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz\n", - " # Generate QC report to assess labeled segmentation\n", - " !sct_qc -i {subject}_acq-CoV_T2starw_crop.nii.gz -s {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -p sct_label_vertebrae -qc {path_qc} -qc-subject {subject}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c46e8d9a", - "metadata": {}, - "outputs": [], - "source": [ - "# Register *_T2starw to CoV_T2starw ⏳\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " for shim_mode in shim_modes:\n", - " # Don't do it for CoV_T2starw\n", - " if shim_mode != 'CoV':\n", - " !sct_register_multimodal -i {subject}_acq-{shim_mode}_T2starw.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -dseg {subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=dl -qc {path_qc}" - ] - }, - { - "cell_type": "markdown", - "id": "aac9d889", - "metadata": {}, - "source": [ - "### Verify QC report (GRE segmentation)\n", - "\n", - "Open the quality control (QC) report located under `ds004906/qc/index.html`. Make sure the spinal cord segmentations are correct before resuming the analysis.\n", - "\n", - ">If you run this notebook on Google Colab, you can skip as the QC report cannot easily be viewed from Google Colab. If you *really* want to see the QC report, you can create a cell where you zip the 'qc/' folder, then you can download it on your local station and open the index.html file. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "90deb5e5", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# Extract the signal intensity on the GRE scan within the spinal cord between levels C3 and T2 (included), which correspond to the region where RF shimming was prescribed\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " for shim_mode in shim_modes:\n", - " # Shim methods are registered to the CoV T2starw scan, so we need to use the added suffix to identify them\n", - " if shim_mode == 'CoV':\n", - " file_suffix = 'crop'\n", - " else:\n", - " file_suffix = 'reg'\n", - " fname_result_sc = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\")\n", - " !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_sc}\n", - " fname_result_csf = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-CSF.csv\")\n", - " !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_csf}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a00909b8", - "metadata": {}, - "outputs": [], - "source": [ - "# Make figure of CSF/SC signal ratio from T2starw scan\n", - "\n", - "# Go back to root data folder\n", - "os.chdir(path_data)\n", - "\n", - "def smooth_data(data, window_size=50):\n", - " \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n", - " return uniform_filter1d(data, size=window_size, mode='nearest')\n", - "\n", - "# Fixed grid for x-axis\n", - "x_grid = np.linspace(0, 1, 100)\n", - "\n", - "# z-slices corresponding to levels C3 to T2 on the PAM50 template. These will be used to scale the x-label of each subject.\n", - "original_vector = np.array([907, 870, 833, 800, 769, 735, 692, 646])\n", - "\n", - "# Normalize the PAM50 z-slice numbers to the 1-0 range (to show inferior-superior instead of superior-inferior)\n", - "min_val = original_vector.min()\n", - "max_val = original_vector.max()\n", - "normalized_vector = 1 - ((original_vector - min_val) / (max_val - min_val))\n", - "\n", - "# Use this normalized vector as x-ticks\n", - "custom_xticks = normalized_vector\n", - "\n", - "# Vertebral level labels\n", - "vertebral_levels = [\"C3\", \"C4\", \"C5\", \"C6\", \"C7\", \"T1\", \"T2\"]\n", - "label_positions = normalized_vector[:-1] + np.diff(normalized_vector) / 2\n", - "\n", - "# Number of subjects determines the number of rows in the subplot\n", - "n_rows = len(subjects)\n", - "\n", - "# Create a figure with multiple subplots\n", - "fig, axes = plt.subplots(n_rows, 1, figsize=(10, 6 * n_rows))\n", - "font_size = 18\n", - "line_width = 3\n", - "\n", - "# Check if axes is an array or a single object\n", - "if n_rows == 1:\n", - " axes = [axes]\n", - "\n", - "# Data storage for statistics\n", - "data_stats = []\n", - "\n", - "# Iterate over each subject and create a subplot\n", - "for i, subject in enumerate(subjects):\n", - " ax = axes[i]\n", - "\n", - " for shim_mode in shim_modes:\n", - " # Initialize list to collect data for this shim method\n", - " method_data = []\n", - "\n", - " # Get signal in SC\n", - " file_csv = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\")\n", - " df = pd.read_csv(file_csv)\n", - " data_sc = df['WA()']\n", - " data_sc_smoothed = smooth_data(data_sc)\n", - "\n", - " # Get signal in CSF\n", - " file_csv = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-CSF.csv\")\n", - " df = pd.read_csv(file_csv)\n", - " data_csf = df['WA()']\n", - " data_csf_smoothed = smooth_data(data_csf)\n", - " \n", - " # Compute ratio\n", - " data_sc_csf_ratio = data_csf_smoothed / data_sc_smoothed\n", - "\n", - " # Normalize the x-axis to a 1-0 scale for each subject (to go from superior-inferior direction)\n", - " x_subject = np.linspace(1, 0, len(data_sc_csf_ratio))\n", - "\n", - " # Interpolate to the fixed grid\n", - " interp_func = interp1d(x_subject, data_sc_csf_ratio, kind='linear', bounds_error=False, fill_value='extrapolate')\n", - " resampled_data = interp_func(x_grid)\n", - "\n", - " method_data.append(resampled_data)\n", - "\n", - " # If there's data for this shim method, plot it and compute stats\n", - " if method_data:\n", - " # Plotting each file's data separately\n", - " for resampled_data in method_data:\n", - " ax.plot(x_grid, resampled_data, label=f\"{shim_mode}\", linewidth=line_width)\n", - " \n", - " # Compute stats on the non-resampled data (to avoid interpolation errors)\n", - " mean_data = np.mean(data_sc_csf_ratio)\n", - " sd_data = np.std(data_sc_csf_ratio)\n", - " data_stats.append([subject, shim_mode, mean_data, sd_data])\n", - "\n", - " # Set x-ticks and labels for the bottom subplot\n", - " if i == n_rows - 1: # Check if it's the last subplot\n", - " # Create a secondary axis for vertebral level labels\n", - " secax = ax.secondary_xaxis('bottom')\n", - " secax.set_xticks(label_positions)\n", - " secax.set_xticklabels(vertebral_levels, fontsize=font_size-4)\n", - " secax.tick_params(axis='x', which='major', length=0) # Hide tick marks\n", - "\n", - " ax.set_xticks(custom_xticks)\n", - " ax.set_xticklabels([''] * len(custom_xticks))\n", - " ax.tick_params(axis='x', which='major', length=0) # Hide tick marks for the primary axis\n", - "\n", - " ax.set_title(f'{subject}', fontsize=font_size)\n", - " ax.set_ylabel('CSF/Cord T2starw signal ratio', fontsize=font_size)\n", - " ax.tick_params(axis='y', which='major', labelsize=font_size-4)\n", - "\n", - " # Add legend only to the first subplot\n", - " if i == 0:\n", - " ax.legend(fontsize=font_size)\n", - "\n", - " ax.grid(True)\n", - "\n", - "# Adjust the layout so labels and titles do not overlap\n", - "plt.tight_layout()\n", - "plt.savefig(os.path.join(path_results, 'fig_gre_csf-sc_ratio.png'), dpi=300, format='png')\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "adbc3046", - "metadata": {}, - "outputs": [], - "source": [ - "# Perform statistics\n", - "\n", - "# Convert to DataFrame and save to CSV\n", - "df_stats = pd.DataFrame(data_stats, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n", - "df_stats.to_csv(os.path.join(path_results, 'stats_t2starw.csv'), index=False)\n", - "\n", - "# Compute statistics across subjects\n", - "grouped_means = df_stats.groupby('Shim_Mode').agg({'Average': ['mean', 'std'], 'Standard_Deviation': ['mean', 'std']})\n", - "\n", - "# Format mean ± standard deviation\n", - "grouped_means['Average_formatted'] = grouped_means['Average']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Average']['std'].map(\"{:.2f}\".format)\n", - "grouped_means['Standard_Deviation_formatted'] = grouped_means['Standard_Deviation']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Standard_Deviation']['std'].map(\"{:.2f}\".format)\n", - "\n", - "# Drop multi-level index and only keep formatted columns\n", - "grouped_means = grouped_means.drop(columns=['Average', 'Standard_Deviation'])\n", - "grouped_means.columns = ['Average', 'Standard_Deviation'] # Rename columns for clarity\n", - "grouped_means.reset_index().to_csv(os.path.join(path_results, 'stats_t2starw_average_across_subjects.csv'), index=False)\n", - "\n", - "# Perform Repeated Measures ANOVA\n", - "aovrm = AnovaRM(df_stats, 'Average', 'Subject', within=['Shim_Mode'])\n", - "res = aovrm.fit()\n", - "print(res.summary())\n", - "\n", - "# Perform Post Hoc paired t-tests\n", - "\n", - "# Filter out subjects that don't have all shim modes\n", - "shim_modes = df_stats['Shim_Mode'].unique()\n", - "valid_subjects = df_stats.groupby('Subject').filter(lambda x: all(mode in x['Shim_Mode'].values for mode in shim_modes))\n", - "\n", - "# Pairwise comparisons\n", - "pairs = list(itertools.combinations(shim_modes, 2))\n", - "p_values = []\n", - "\n", - "for pair in pairs:\n", - " data1 = valid_subjects[valid_subjects['Shim_Mode'] == pair[0]]['Average']\n", - " data2 = valid_subjects[valid_subjects['Shim_Mode'] == pair[1]]['Average']\n", - " _, p_value = ttest_rel(data1, data2)\n", - " p_values.append(p_value)\n", - "\n", - "# Function for Benjamini-Hochberg FDR correction\n", - "def benjamini_hochberg(p_values):\n", - " n = len(p_values)\n", - " sorted_p_values = np.sort(p_values)\n", - " sorted_index = np.argsort(p_values)\n", - " adjusted_p_values = np.zeros(n)\n", - "\n", - " for i, p in enumerate(sorted_p_values):\n", - " adjusted_p_values[sorted_index[i]] = min(p * n / (i + 1), 1)\n", - "\n", - " return adjusted_p_values\n", - "\n", - "# Applying Benjamini-Hochberg FDR correction\n", - "adjusted_p_values = benjamini_hochberg(p_values)\n", - "\n", - "# Creating a summary DataFrame\n", - "comparison_results = pd.DataFrame({'Pair': pairs, 'P-Value': p_values, 'Adjusted P-Value': adjusted_p_values})\n", - "\n", - "# Save to CSV\n", - "comparison_results.to_csv(os.path.join(path_results, 'stats_t2starw_paired_ttest.csv'), index=False)\n", - "\n", - "print(f\"Paired t-tests:\\n{comparison_results}\")" - ] - }, - { - "cell_type": "markdown", - "id": "bda7c3ac", - "metadata": {}, - "source": [ - "## Process fmap/TFL (flip angle maps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e0cde907", - "metadata": {}, - "outputs": [], - "source": [ - "# Register TFL flip angle maps to the GRE scan ⏳\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " !sct_register_multimodal -i {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -d ../anat/{subject}_acq-CoV_T2starw_crop.nii.gz -dseg ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc {path_qc}" - ] - }, - { - "cell_type": "markdown", - "id": "39de4623", - "metadata": {}, - "source": [ - "### Verify QC report (B1maps to GRE registration)\n", - "\n", - "Open the QC report located under `ds004906/qc/index.html`. Make sure the registration are correct before resuming the analysis." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6a6d7670", - "metadata": {}, - "outputs": [], - "source": [ - "# Warping spinal cord segmentation and vertebral level to each flip angle map\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " # Use linear interpolation to preserve partial volume information (needed when extracting values along the cord)\n", - " !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x linear -o {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz\n", - " # Use nearest neighbour (nn) interpolation because we are dealing with non-binary discrete labels\n", - " !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x nn -o {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a835cdee", - "metadata": {}, - "outputs": [], - "source": [ - "# Convert the flip angle maps to B1+ efficiency maps [nT/V] (inspired by code from Kyle Gilbert)\n", - "# The approach consists in calculating the B1+ efficiency using a 1ms, pi-pulse at the acquisition voltage,\n", - "# then scale the efficiency by the ratio of the measured flip angle to the requested flip angle in the pulse sequence.\n", - "\n", - "GAMMA = 2.675e8; # [rad / (s T)]\n", - "requested_fa = 90 # saturation flip angle -- hard-coded in sequence\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " # Fetch the reference voltage from the JSON sidecar to the TFL B1map sequence\n", - " with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n", - " metadata = json.load(f)\n", - " ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n", - " print(f\"ref_voltage [V]: {ref_voltage} ({subject}_acq-famp{shim_mode}_TB1TFL)\")\n", - "\n", - " # Open flip angle map with nibabel\n", - " nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n", - " acquired_fa = nii.get_fdata()\n", - "\n", - " # Siemens maps are in units of flip angle * 10 (in degrees)\n", - " acquired_fa = acquired_fa / 10\n", - "\n", - " # Account for the power loss between the coil and the socket. That number was given by Siemens.\n", - " voltage_at_socket = ref_voltage * 10 ** -0.095\n", - "\n", - " # Compute B1 map in [T/V]\n", - " b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n", - "\n", - " # Convert to [nT/V]\n", - " b1_map = b1_map * 1e9\n", - "\n", - " # Save as NIfTI file\n", - " nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n", - " nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "740fda91", - "metadata": {}, - "outputs": [], - "source": [ - "# Extract B1+ value along the spinal cord between levels C3 and T2 (included)\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " fname_result_b1plus = os.path.join(path_results, f\"{subject}_TB1map_{shim_mode}.csv\")\n", - " !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o {fname_result_b1plus}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b0185bb0", - "metadata": {}, - "outputs": [], - "source": [ - "# Make figure of B1+ values along the spinal cord across shim methods\n", - "\n", - "# Go back to root data folder\n", - "os.chdir(os.path.join(path_data))\n", - "\n", - "def smooth_data(data, window_size=20):\n", - " \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n", - " return uniform_filter1d(data, size=window_size, mode='nearest')\n", - "\n", - "# Fixed grid for x-axis\n", - "x_grid = np.linspace(0, 1, 100)\n", - "\n", - "# z-slices corresponding to levels C3 to T2 on the PAM50 template. These will be used to scale the x-label of each subject.\n", - "original_vector = np.array([907, 870, 833, 800, 769, 735, 692, 646])\n", - "\n", - "# Normalize the PAM50 z-slice numbers to the 1-0 range (to show inferior-superior instead of superior-inferior)\n", - "min_val = original_vector.min()\n", - "max_val = original_vector.max()\n", - "normalized_vector = 1 - ((original_vector - min_val) / (max_val - min_val))\n", - "\n", - "# Use this normalized vector as x-ticks\n", - "custom_xticks = normalized_vector\n", - "\n", - "# Vertebral level labels\n", - "vertebral_levels = [\"C3\", \"C4\", \"C5\", \"C6\", \"C7\", \"T1\", \"T2\"]\n", - "# Calculate midpoints for label positions\n", - "label_positions = normalized_vector[:-1] + np.diff(normalized_vector) / 2\n", - "\n", - "# Number of subjects determines the number of rows in the subplot\n", - "n_rows = len(subjects)\n", - "\n", - "# Create a figure with multiple subplots\n", - "fig, axes = plt.subplots(n_rows, 1, figsize=(10, 6 * n_rows))\n", - "font_size = 18\n", - "line_width = 3\n", - "\n", - "# Check if axes is an array or a single object\n", - "if n_rows == 1:\n", - " axes = [axes]\n", - "\n", - "# Data storage for statistics\n", - "data_stats = []\n", - "\n", - "# Iterate over each subject and create a subplot\n", - "for i, subject in enumerate(subjects):\n", - " ax = axes[i]\n", - " \n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - "\n", - " for shim_mode in shim_modes:\n", - " # Initialize list to collect data for this shim method\n", - " method_data = []\n", - "\n", - " file_csv = os.path.join(path_results, f\"{subject}_TB1map_{shim_mode}.csv\")\n", - " df = pd.read_csv(file_csv)\n", - " wa_data = df['WA()']\n", - "\n", - " # Normalize the x-axis to a 1-0 scale for each subject (to go from superior-inferior direction)\n", - " x_subject = np.linspace(1, 0, len(wa_data))\n", - "\n", - " # Interpolate to the fixed grid\n", - " interp_func = interp1d(x_subject, wa_data, kind='linear', bounds_error=False, fill_value='extrapolate')\n", - " resampled_data = interp_func(x_grid)\n", - "\n", - " # Apply smoothing\n", - " smoothed_data = smooth_data(resampled_data)\n", - "\n", - " method_data.append(smoothed_data)\n", - "\n", - " # If there's data for this shim method, plot it\n", - " if method_data:\n", - " # Plotting each file's data separately\n", - " for resampled_data in method_data:\n", - " ax.plot(x_grid, resampled_data, label=f\"{shim_mode}\", linewidth=line_width)\n", - "\n", - " # Compute stats on the non-resampled data (to avoid interpolation errors)\n", - " mean_data = np.mean(wa_data)\n", - " sd_data = np.std(wa_data)\n", - " data_stats.append([subject, shim_mode, mean_data, sd_data])\n", - "\n", - " # Set x-ticks and labels for the bottom subplot\n", - " if i == n_rows - 1: # Check if it's the last subplot\n", - " # Create a secondary axis for vertebral level labels\n", - " secax = ax.secondary_xaxis('bottom')\n", - " secax.set_xticks(label_positions)\n", - " secax.set_xticklabels(vertebral_levels, fontsize=font_size-4)\n", - " secax.tick_params(axis='x', which='major', length=0) # Hide tick marks\n", - "\n", - " ax.set_xticks(custom_xticks)\n", - " ax.set_xticklabels([''] * len(custom_xticks))\n", - " ax.tick_params(axis='x', which='major', length=0) # Hide tick marks for the primary axis\n", - "\n", - " ax.set_title(f'{subject}', fontsize=font_size)\n", - " ax.set_ylabel('B1+ efficiency [nT/V]', fontsize=font_size)\n", - " ax.tick_params(axis='y', which='major', labelsize=font_size-4)\n", - "\n", - " # Add legend only to the first subplot\n", - " if i == 0:\n", - " ax.legend(fontsize=font_size)\n", - "\n", - " ax.grid(True)\n", - "\n", - "# Adjust the layout so labels and titles do not overlap\n", - "plt.tight_layout()\n", - "plt.savefig(os.path.join(path_results, 'fig_b1plus.png'), dpi=300, format='png')\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c09f3c05", - "metadata": {}, - "outputs": [], - "source": [ - "# Perform statistics\n", - "\n", - "# Convert to DataFrame and save to CSV\n", - "df_stats = pd.DataFrame(data_stats, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n", - "df_stats.to_csv(os.path.join(path_results, 'stats_b1plus.csv'), index=False)\n", - "\n", - "# Compute statistics across subjects\n", - "grouped_means = df_stats.groupby('Shim_Mode').agg({'Average': ['mean', 'std'], 'Standard_Deviation': ['mean', 'std']})\n", - "\n", - "# Format mean ± standard deviation\n", - "grouped_means['Average_formatted'] = grouped_means['Average']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Average']['std'].map(\"{:.2f}\".format)\n", - "grouped_means['Standard_Deviation_formatted'] = grouped_means['Standard_Deviation']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Standard_Deviation']['std'].map(\"{:.2f}\".format)\n", - "\n", - "# Drop multi-level index and only keep formatted columns\n", - "grouped_means = grouped_means.drop(columns=['Average', 'Standard_Deviation'])\n", - "grouped_means.columns = ['Average', 'Standard_Deviation'] # Rename columns for clarity\n", - "grouped_means.reset_index().to_csv(os.path.join(path_results, 'stats_b1plus_average_across_subjects.csv'), index=False)\n", - "\n", - "# Perform Repeated Measures ANOVA\n", - "aovrm = AnovaRM(df_stats, 'Average', 'Subject', within=['Shim_Mode'])\n", - "res = aovrm.fit()\n", - "print(res.summary())\n", - "\n", - "# Perform Post Hoc paired t-tests\n", - "\n", - "# Filter out subjects that don't have all shim modes\n", - "shim_modes = df_stats['Shim_Mode'].unique()\n", - "valid_subjects = df_stats.groupby('Subject').filter(lambda x: all(mode in x['Shim_Mode'].values for mode in shim_modes))\n", - "\n", - "# Pairwise comparisons\n", - "pairs = list(itertools.combinations(shim_modes, 2))\n", - "p_values = []\n", - "\n", - "for pair in pairs:\n", - " data1 = valid_subjects[valid_subjects['Shim_Mode'] == pair[0]]['Average']\n", - " data2 = valid_subjects[valid_subjects['Shim_Mode'] == pair[1]]['Average']\n", - " _, p_value = ttest_rel(data1, data2)\n", - " p_values.append(p_value)\n", - "\n", - "# Function for Benjamini-Hochberg FDR correction\n", - "def benjamini_hochberg(p_values):\n", - " n = len(p_values)\n", - " sorted_p_values = np.sort(p_values)\n", - " sorted_index = np.argsort(p_values)\n", - " adjusted_p_values = np.zeros(n)\n", - "\n", - " for i, p in enumerate(sorted_p_values):\n", - " adjusted_p_values[sorted_index[i]] = min(p * n / (i + 1), 1)\n", - "\n", - " return adjusted_p_values\n", - "\n", - "# Applying Benjamini-Hochberg FDR correction\n", - "adjusted_p_values = benjamini_hochberg(p_values)\n", - "\n", - "# Creating a summary DataFrame\n", - "comparison_results = pd.DataFrame({'Pair': pairs, 'P-Value': p_values, 'Adjusted P-Value': adjusted_p_values})\n", - "\n", - "# Save to CSV\n", - "comparison_results.to_csv(os.path.join(path_results, 'stats_b1plus_paired_ttest.csv'), index=False)\n", - "\n", - "print(f\"Paired t-tests:\\n{comparison_results}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1dde4c00", - "metadata": {}, - "outputs": [], - "source": [ - "# Prepare RF shimming mask for figure\n", - "\n", - "# Select subject to show\n", - "subject = 'sub-05'\n", - "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - "\n", - "# Create RF shimming mask from the segmentation (re-create the live RF shimming procedure)\n", - "file_mask = f\"{subject}_acq-anatCP_TB1TFL_mask-shimming.nii.gz\"\n", - "!sct_maths -i {subject}_acq-anatCP_TB1TFL_seg_labeled.nii.gz -thr 3 -uthr 9 -o {file_mask}\n", - "!sct_maths -i {file_mask} -bin 1 -o {file_mask}\n", - "!sct_create_mask -i {subject}_acq-anatCP_TB1TFL.nii.gz -p centerline,{file_mask} -size 28mm -f cylinder -o {file_mask}\n", - "\n", - "# Dilate mask\n", - "!sct_maths -i {file_mask} -dilate 2 -shape disk -dim 2 -o {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz\n", - "# Subtract dilated mask with mask to get the edge\n", - "!sct_maths -i {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz -sub {file_mask} -o {file_mask}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "93d72ec3", - "metadata": {}, - "outputs": [], - "source": [ - "# Create figure of B1+ maps\n", - "\n", - "# Select subject to show\n", - "subject = 'sub-05'\n", - "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - "file_mask = f\"{subject}_acq-anatCP_TB1TFL_mask-shimming.nii.gz\"\n", - "\n", - "# Load the statistics from the CSV file\n", - "stats_df = pd.read_csv(os.path.join(path_results, 'stats_b1plus.csv'))\n", - "\n", - "# Filter for the specific subject\n", - "subject_stats = stats_df[stats_df['Subject'] == subject]\n", - "\n", - "# Defining crop limits for resulting figure\n", - "xmin = 20\n", - "xmax = 110\n", - "ymin = 20\n", - "ymax = 150\n", - "\n", - "# Defining dynamic range\n", - "dynmin = 5 \n", - "dynmax = 30\n", - "\n", - "# Create a figure with multiple subplots\n", - "fig, axes = plt.subplots(2, 4, figsize=(8, 6))\n", - "font_size = 14\n", - "axes=axes.flatten()\n", - "\n", - "# First, plot the anatomical image with an overlay of the mask\n", - "\n", - "# Load data\n", - "CP_anat=nib.load(f\"{subject}_acq-anatCP_TB1TFL.nii.gz\")\n", - "CP_SC=nib.load(file_mask)\n", - "CP_nTpV=nib.load(f\"{subject}_acq-CP_TB1map.nii.gz\")\n", - "\n", - "# Defining mask based on the magnitude image intensity threshold\n", - "cslice=CP_anat.shape[2] // 2 -2 #shows the SC seg best\n", - "threshold=300\n", - "mask=CP_anat.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n", - "mask=np.where(mask > threshold, 1, 0)\n", - "\n", - "# Cropping anat, SC, B1+\n", - "CP_anat=CP_anat.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n", - "CP_anat=CP_anat*mask\n", - "CP_SC=CP_SC.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n", - "CP_nTpV=CP_nTpV.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n", - "CP_nTpV=CP_nTpV*mask\n", - "\n", - "# All opacity overalys look ugly: workaround, set the anat slice to a max value where the segmentation exists\n", - "CP_anat[CP_SC>0.5]=4000;\n", - "\n", - "# Plotting anat overlayed with SC\n", - "splot=axes[0]\n", - "splot.imshow((CP_anat.T), cmap='gray', origin='lower',vmin=0,vmax=2000)#, interpolation='spline36')\n", - "splot.set_title('Anat', size=font_size)\n", - "splot.axis('off')\n", - "\n", - "# Then, plot each B1+ map, with an overlay of the mean and CV inside the cord\n", - "for i,shim_mode in enumerate(shim_modes):\n", - " # Load data\n", - " B1map=nib.load(f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")\n", - " B1map=B1map.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n", - " B1map=B1map*mask\n", - "\n", - " # Plot\n", - " splot=axes[i+1]\n", - " im = splot.imshow((B1map.T), cmap='viridis', origin='lower',vmin=dynmin,vmax=dynmax)#, interpolation='spline36')\n", - " splot.set_title(shim_mode, size=font_size)\n", - " splot.axis('off')\n", - "\n", - " # Find the statistics for the current shim mode\n", - " shim_stats = subject_stats[subject_stats['Shim_Mode'] == shim_mode]\n", - " if not shim_stats.empty:\n", - " mean_val = shim_stats.iloc[0]['Average']\n", - " std_val = shim_stats.iloc[0]['Standard_Deviation']\n", - " cv = std_val / mean_val * 100 # Coefficient of variation in percentage\n", - " annotation_text = f\"{mean_val:.2f} nT/V\\n{cv:.2f}%\"\n", - " splot.annotate(annotation_text, (0.05, 0.95), xycoords='axes fraction', \n", - " fontsize=10, color='white', \n", - " verticalalignment='top', horizontalalignment='left')\n", - "\n", - "plt.tight_layout()\n", - "plt.subplots_adjust(wspace=0.1, hspace=0, right=0.9)\n", - "\n", - "# Colorbar\n", - "# Assume that the colorbar should start at the bottom of the lower row of subplots and\n", - "# extend to the top of the upper row of subplots\n", - "cbar_bottom = 0.06 # This might need adjustment\n", - "cbar_height = 0.83 # This represents the total height of both rows of subplots\n", - "cbar_ax = fig.add_axes([0.93, cbar_bottom, 0.03, cbar_height])\n", - "cbar = plt.colorbar(im, cax=cbar_ax)\n", - "\n", - "# cbar_ax = fig.add_axes([0.95, 0.5, 0.04, 0.4])\n", - "# cbar = plt.colorbar(im, cax=cbar_ax)\n", - "cbar_ax.set_title('nT/V', size=12)\n", - "plt.savefig(os.path.join(path_results, 'fig_b1plus_map.png'), dpi=300, format='png')\n", - "plt.show\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9128d799", - "metadata": {}, - "outputs": [], - "source": [ - "# Indicate duration of data processing\n", - "\n", - "end_time = datetime.now()\n", - "total_time = (end_time - start_time).total_seconds()\n", - "\n", - "# Convert seconds to a timedelta object\n", - "total_time_delta = timedelta(seconds=total_time)\n", - "\n", - "# Format the timedelta object to a string\n", - "formatted_time = str(total_time_delta)\n", - "\n", - "# Pad the string representation if less than an hour\n", - "formatted_time = formatted_time.rjust(8, '0')\n", - "\n", - "print(f\"Total Runtime [hour:min:sec]: {formatted_time}\")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.8" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -}