Volume-wise T2*/S0 estimation with t2smap

Volume-wise T2*/S0 estimation with t2smap#

Use tedana.workflows.t2smap_workflow() [DuPre et al., 2021] to calculate volume-wise T2*/S0, as in Power et al. [2018] and Heunis et al. [2021].

import json
import os
from glob import glob

import matplotlib.pyplot as plt
import nibabel as nb
import numpy as np
from IPython import display
from myst_nb import glue
from nilearn import image, plotting
from tedana import workflows

data_path = os.path.abspath('../DATA')
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
func_dir = os.path.join(data_path, "ds006185/sub-24053/ses-1/func/")
data_files = sorted(
    glob(
        os.path.join(
            func_dir,
            "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_echo-*_part-mag_desc-preproc_bold.nii.gz",
        ),
    ),
)
echo_times = []
for f in data_files:
    json_file = f.replace('.nii.gz', '.json')
    with open(json_file, 'r') as fo:
        metadata = json.load(fo)
    echo_times.append(metadata['EchoTime'] * 1000)
mask_file = os.path.join(
    func_dir,
    "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_part-mag_desc-brain_mask.nii.gz"
)
confounds_file = os.path.join(
    func_dir,
    "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_part-mag_desc-confounds_timeseries.tsv",
)

out_dir = os.path.join(data_path, "fit")
workflows.t2smap_workflow(
    data_files,
    echo_times,
    out_dir=out_dir,
    mask=mask_file,
    prefix="sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_",
    fittype="loglin",
    fitmode="ts",
    overwrite=True,
)
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[3], line 1
----> 1 workflows.t2smap_workflow(
      2     data_files,
      3     echo_times,
      4     out_dir=out_dir,
      5     mask=mask_file,
      6     prefix="sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_",
      7     fittype="loglin",
      8     fitmode="ts",
      9     overwrite=True,
     10 )

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/tedana/workflows/t2smap.py:322, in t2smap_workflow(data, tes, n_independent_echos, out_dir, mask, prefix, convention, dummy_scans, exclude, masktype, fittype, fitmode, combmode, debug, verbose, quiet, overwrite, t2smap_command)
    320 out_dir = op.abspath(out_dir)
    321 if not op.isdir(out_dir):
--> 322     os.mkdir(out_dir)
    324 # Parse exclude parameter
    325 exclude_idx = parse_volume_indices(exclude)

FileNotFoundError: [Errno 2] No such file or directory: '/home/runner/work/multi-echo-data-analysis/multi-echo-data-analysis/DATA/fit'
out_files = sorted(glob(os.path.join(out_dir, "*")))
out_files = [os.path.basename(f) for f in out_files]
print("\n".join(out_files))

Hide code cell content

fig, axes = plt.subplots(figsize=(16, 16), nrows=3)

plotting.plot_carpet(
    os.path.join(out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_desc-optcom_bold.nii.gz"),
    axes=axes[0],
    figure=fig,
)
axes[0].set_title("Optimally Combined Data", fontsize=20)
plotting.plot_carpet(
    os.path.join(out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_T2starmap.nii.gz"),
    axes=axes[1],
    figure=fig,
)
axes[1].set_title("T2* Estimates", fontsize=20)
plotting.plot_carpet(
    os.path.join(out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_S0map.nii.gz"),
    axes=axes[2],
    figure=fig,
)
axes[2].set_title("S0 Estimates", fontsize=20)
axes[0].xaxis.set_visible(False)
axes[1].xaxis.set_visible(False)
axes[0].spines["bottom"].set_visible(False)
axes[1].spines["bottom"].set_visible(False)
fig.tight_layout()

glue("figure_volumewise_t2ss0_carpets", fig, display=False)
fig, ax = plt.subplots(figsize=(16, 8))
plotting.plot_stat_map(
    image.mean_img(
        os.path.join(out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_T2starmap.nii.gz")
    ),
    vmax=0.6,
    draw_cross=False,
    bg_img=None,
    figure=fig,
    axes=ax,
)
glue("figure_mean_volumewise_t2s", fig, display=True)
fig, ax = plt.subplots(figsize=(16, 8))
plotting.plot_stat_map(
    image.mean_img(
        os.path.join(out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_S0map.nii.gz")
    ),
    vmax=8000,
    draw_cross=False,
    bg_img=None,
    figure=fig,
    axes=ax,
)
glue("figure_mean_volumewise_s0", fig, display=True)
fig, axes = plt.subplots(figsize=(16, 15), nrows=5)
plotting.plot_epi(
    image.mean_img(data_files[0]),
    draw_cross=False,
    bg_img=None,
    cut_coords=[-10, 0, 10, 20, 30, 40, 50, 60, 70],
    display_mode="z",
    axes=axes[0],
)
plotting.plot_epi(
    image.mean_img(data_files[1]),
    draw_cross=False,
    bg_img=None,
    cut_coords=[-10, 0, 10, 20, 30, 40, 50, 60, 70],
    display_mode="z",
    axes=axes[1],
)
plotting.plot_epi(
    image.mean_img(data_files[2]),
    draw_cross=False,
    bg_img=None,
    cut_coords=[-10, 0, 10, 20, 30, 40, 50, 60, 70],
    display_mode="z",
    axes=axes[2],
)
plotting.plot_epi(
    image.mean_img(data_files[3]),
    draw_cross=False,
    bg_img=None,
    cut_coords=[-10, 0, 10, 20, 30, 40, 50, 60, 70],
    display_mode="z",
    axes=axes[3],
)
plotting.plot_epi(
    image.mean_img(
        os.path.join(
            out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_desc-optcom_bold.nii.gz"
        )
    ),
    draw_cross=False,
    bg_img=None,
    cut_coords=[-10, 0, 10, 20, 30, 40, 50, 60, 70],
    display_mode="z",
    axes=axes[4],
)
glue("figure_mean_echos_and_optcom", fig, display=True)
te30_tsnr = image.math_img(
    "(np.nanmean(img, axis=3) / np.nanstd(img, axis=3)) * mask",
    img=data_files[1],
    mask=mask_file,
)
oc_tsnr = image.math_img(
    "(np.nanmean(img, axis=3) / np.nanstd(img, axis=3)) * mask",
    img=os.path.join(
        out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_desc-optcom_bold.nii.gz"
    ),
    mask=mask_file,
)
vmax = np.nanmax(np.abs(oc_tsnr.get_fdata()))

fig, axes = plt.subplots(figsize=(10, 8), nrows=2)
plotting.plot_stat_map(
    te30_tsnr,
    draw_cross=False,
    bg_img=None,
    threshold=0.1,
    cut_coords=[0, 10, 10],
    vmax=vmax,
    symmetric_cbar=False,
    axes=axes[0],
)
axes[0].set_title("TE30 TSNR", fontsize=16)
plotting.plot_stat_map(
    oc_tsnr,
    draw_cross=False,
    bg_img=None,
    threshold=0.1,
    cut_coords=[0, 10, 10],
    vmax=vmax,
    symmetric_cbar=False,
    axes=axes[1],
)
axes[1].set_title("Optimal Combination TSNR", fontsize=16)
glue("figure_t2snr_te30_and_optcom", fig, display=True)
fig, ax = plt.subplots(figsize=(16, 8))
plotting.plot_carpet(
    data_files[1],
    axes=ax,
)
glue("figure_echo3_carpet", fig, display=True)
fig, ax = plt.subplots(figsize=(16, 8))
plotting.plot_carpet(
    os.path.join(out_dir, "sub-24053_ses-1_task-rat_rec-nordic_dir-PA_run-01_desc-optcom_bold.nii.gz"),
    axes=ax,
)
glue("figure_carpet_optcom", fig, display=True)