NiWrap and Connectome Workbench

# Load connectomeworkbench
import module
await module.load('connectomeworkbench/2.1.0')
await module.list()

%%capture 
!pip install niwrap==0.6.1 nilearn==0.12.0 nibabel==5.3.2 pandas==2.2.3

# Import the necessary libraries

import matplotlib.pyplot as plt 
import os 
from pathlib import Path
from niwrap_workbench.workbench import (cifti_math, cifti_math_var_params, cifti_separate, cifti_separate_metric_params, 
                                        cifti_stats, metric_smoothing, cifti_create_dense_scalar, cifti_create_dense_scalar_left_metric_params, cifti_create_dense_scalar_right_metric_params,
                                        cifti_parcellate, file_information, cifti_label_export_table)
import re
import pandas as pd
from nilearn import plotting, surface
import nibabel as nib
import numpy as np

%%bash
# --- Configuration ---
REPO_URL="https://github.com/OpenNeuroDatasets/ds004496.git"
DATASET_DIR="ds004496"
SUBJECT_ID="sub-01"

# --- Datalad Installation ---
if [ ! -d "$DATASET_DIR" ]; then
  datalad install $REPO_URL
else
  echo "Dataset directory '$DATASET_DIR' already exists. Skipping install."
fi

cd $DATASET_DIR

# --- Datalad Data Download ---
echo "Downloading Ciftify derivative files for $SUBJECT_ID from ds004496..."

# 1. Get the beta map from the floc task.
datalad get "derivatives/ciftify/$SUBJECT_ID/results/ses-floc_task-floc/ses-floc_task-floc_beta.dscalar.nii"

# 2. Get the corresponding T1 surface geometry files.
#datalad get "derivatives/ciftify/$SUBJECT_ID/T1w_fsLR_surface"

# 2. Get the corresponding standard surface geometry files.
datalad get "derivatives/ciftify/$SUBJECT_ID/standard_fsLR_surface"

echo "Datalad download complete for ds004496."

# --- Define File Paths ---
DATASET_DIR = Path("ds004496")
SUBJECT_ID = "sub-01"

# The root directory of the Ciftify derivatives
deriv_dir = DATASET_DIR / "derivatives" / "ciftify"

# 1. Input File: The statistical map
beta_map_path = (
    deriv_dir
    / SUBJECT_ID
    / "results"
    / "ses-floc_task-floc"
    / "ses-floc_task-floc_beta.dscalar.nii"
)

# 2. Input Files: The surface geometry files.
surface_dir = deriv_dir / SUBJECT_ID / "standard_fsLR_surface"
left_surf_path = surface_dir / f"{SUBJECT_ID}.L.midthickness.32k_fs_LR.surf.gii"
right_surf_path = surface_dir / f"{SUBJECT_ID}.R.midthickness.32k_fs_LR.surf.gii"
left_infl_path = surface_dir / f"{SUBJECT_ID}.L.inflated.32k_fs_LR.surf.gii"
right_infl_path = surface_dir / f"{SUBJECT_ID}.R.inflated.32k_fs_LR.surf.gii"
sulcal_depth_path = surface_dir / f"{SUBJECT_ID}.sulc.32k_fs_LR.dscalar.nii"


# --- Verify files exist before we start ---
assert left_surf_path.exists(), f"Left mid-surface not found: {left_surf_path}"
assert right_surf_path.exists(), f"Right mid-surface not found: {right_surf_path}"
assert left_infl_path.exists(), f"Left inflated surface not found: {left_surf_path}"
assert right_infl_path.exists(), f"Right inflated surface not found: {right_surf_path}"
assert sulcal_depth_path.exists(), f"Sulcal depth not found: {sulcal_depth_path}"
print("✅ All necessary files found. Ready for analysis.")

# Download Schaefer atlas (skip if exists)
! wget -nc -q --show-progress https://github.com/ThomasYeoLab/CBIG/raw/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/Parcellations/HCP/fslr32k/cifti/Schaefer2018_400Parcels_17Networks_order.dlabel.nii

#create output directory
output_dir = Path("./niwrap_results").absolute()
output_dir.mkdir(parents=True, exist_ok=True)

! wb_command -file-information $beta_map_path -only-map-names

! wb_command -cifti-merge ./niwrap_results/beta_face.dscalar.nii -cifti $beta_map_path -column 3

! wb_command -file-information ./niwrap_results/beta_face.dscalar.nii -only-map-names

# Load the beta map
img = nib.load("./niwrap_results/beta_face.dscalar.nii")
data = img.get_fdata()

# Mask out NaNs (some CIFTI data may have them)
masked_data = data[~np.isnan(data)]

# Compute 95th percentile
percentile_95 = np.percentile(masked_data, 95)

print(f"95th percentile value: {percentile_95}")

output_file = output_dir / "beta_face_thresholded.dscalar.nii"

result_cifti = cifti_math(
    expression=f"x * (x > {percentile_95})",  # Retain original values > 95th percentile
    cifti_out=str(output_file),
    var=[
        cifti_math_var_params(
            name="x",
            cifti=str("./niwrap_results/beta_face.dscalar.nii")
        )
    ]
)

metric = [
    cifti_separate_metric_params(
        structure="CORTEX_LEFT",
        metric_out=str(output_dir / "sub-01_thresh_lh.func.gii"),
    ),
    cifti_separate_metric_params(
        structure="CORTEX_RIGHT",
        metric_out=str(output_dir / "sub-01_thresh_rh.func.gii"),
    )
]

result_separate = cifti_separate(
    cifti_in=str(result_cifti.cifti_out), #thresholded image
    direction="COLUMN",
    metric=metric
)

metric = [
    cifti_separate_metric_params(
        structure="CORTEX_LEFT",
        metric_out=str(output_dir / "sub-01.L.sulc.func.gii"),
    ),
    cifti_separate_metric_params(
        structure="CORTEX_RIGHT",
        metric_out=str(output_dir / "sub-01.R.sulc.func.gii"),
    )
]

result_separate_sulc = cifti_separate(
    cifti_in=str(sulcal_depth_path),
    direction="COLUMN",
    metric=metric
)

# Plot thresholde maps on cortical surfaces with interactive rotation and zooming

# Load meshes and concatenate coordinates
mesh_lh = surface.load_surf_mesh(str(left_surf_path))
mesh_rh = surface.load_surf_mesh(str(right_surf_path))
coords_combined = np.vstack([mesh_lh.coordinates, mesh_rh.coordinates])

# Offset RH face indices by LH vertex count
offset = mesh_lh.coordinates.shape[0]
faces_rh_offset = mesh_rh.faces + offset

# Concatenate faces
faces_combined = np.vstack([mesh_lh.faces, faces_rh_offset])

# Combine into mesh tuple
combined_mesh_mid = (coords_combined, faces_combined)

# Load and concatenate data
data_lh = nib.load(str(result_separate.metric[0].metric_out)).darrays[0].data
data_rh = nib.load(str(result_separate.metric[1].metric_out)).darrays[0].data
data_combined = np.concatenate([data_lh, data_rh])

# Plot both hemispheres
view = plotting.view_surf(
    surf_mesh=combined_mesh_mid,
    surf_map=data_combined,
    hemi='both',
    threshold=0.0001, #for better contrast of the brain surface
    title="Face > Others: fLoc Functional Localizer Results (Top 5% Activation)",
    title_fontsize=18,
    colorbar=False,
    darkness=None,   
)

view  

#inflated surface
infl_mesh_lh = surface.load_surf_mesh(str(left_infl_path))
infl_mesh_rh = surface.load_surf_mesh(str(right_infl_path))
coords_infl_combined = np.vstack([infl_mesh_lh.coordinates, infl_mesh_rh.coordinates])

facesinfl_rh_offset = infl_mesh_rh.faces + infl_mesh_lh.coordinates.shape[0]
faces_infl_combined = np.vstack([infl_mesh_lh.faces, facesinfl_rh_offset])
combined_infl_mesh = (coords_infl_combined, faces_infl_combined)

#sulcal depth background for anatomical context
sulc_lh = nib.load(str(result_separate_sulc.metric[0].metric_out)).darrays[0].data
sulc_rh = nib.load(str(result_separate_sulc.metric[1].metric_out)).darrays[0].data
sulc_combined = np.concatenate([sulc_lh, sulc_rh])


fig = plotting.plot_surf_stat_map(
    surf_mesh=combined_infl_mesh,
    stat_map=data_combined,
    hemi='both',
    title="Face > Others: fLoc Results (Top 5% Activation) from ventral",
    colorbar=True,
    darkness=None,
    threshold=0.0001,
    bg_map=sulc_combined,
    view='ventral',
    cmap='hot'
)
fig.show()

result_stats = cifti_stats(
    cifti_in=str(result_cifti.cifti_out),
    opt_reduce_operation="COUNT_NONZERO",
    opt_show_map_name=True
)

# Left hemisphere
smoothed_lh = output_dir / "sub-01_thresh_smooth.L.func.gii"

smooth_left = metric_smoothing(
    surface=str(left_surf_path),
    metric_in=str(result_separate.metric[0].metric_out),
    smoothing_kernel=2.0,
    metric_out=str(smoothed_lh),
)

# Right hemisphere
smoothed_rh = output_dir / "sub-01_thresh_smooth.R.func.gii"

smooth_right = metric_smoothing(
    surface=str(right_surf_path),
    metric_in=str(result_separate.metric[1].metric_out),
    smoothing_kernel=2.0,
    metric_out=str(smoothed_rh),
)

# Load and concatenate smoothed data
data_lh_smooth = nib.load(str(smoothed_lh)).darrays[0].data
data_rh_smooth = nib.load(str(smoothed_rh)).darrays[0].data
data_smooth = np.concatenate([data_lh_smooth, data_rh_smooth])


fig = plotting.plot_surf_stat_map(
    surf_mesh=combined_infl_mesh,
    stat_map=data_smooth,
    hemi='both',
    title="Face > Others: fLoc smoothed Results (Top 5% Activation) from ventral",
    colorbar=True,
    darkness=None,
    threshold=0.0001,
    bg_map=sulc_combined,
    view='ventral',
    cmap='hot',
)
fig.show()

parcellated = output_dir / "sub-01_parcellated.pscalar.nii"

result_parc = cifti_parcellate(
    cifti_in=result_cifti.cifti_out,
    cifti_label="Schaefer2018_400Parcels_17Networks_order.dlabel.nii",
    direction="COLUMN",
    cifti_out=str(parcellated)
)

export_table = cifti_label_export_table(
	label_in= str("./Schaefer2018_400Parcels_17Networks_order.dlabel.nii"),
    map_="parcels", 
	table_out= str(output_dir / "schaefer_labels.txt")
)

label_file = "./niwrap_results/schaefer_labels.txt"

with open(label_file, 'r') as f:
    lines = [line.strip() for line in f if line.strip()]

# Separate odd and even lines
names = lines[0::2]   # 0, 2, 4 ... parcel names
colors = lines[1::2]  # 1, 3, 5 ... indices and RGB info

# Parse colors info lines into components
indices = []
Rs, Gs, Bs, As = [], [], [], []

for c in colors:
    # Split by whitespace
    parts = re.split(r'\s+', c)
    # First part is index, next are RGBA
    indices.append(int(parts[0]))
    Rs.append(float(parts[1]))
    Gs.append(float(parts[2]))
    Bs.append(float(parts[3]))
    As.append(float(parts[4]) if len(parts) > 4 else None)



df_labels = pd.DataFrame({
    'Index': indices,
    'Name': names,
    'R': Rs,
    'G': Gs,
    'B': Bs,
    'A': As
})

df_labels.head()

img = nib.load(result_parc.cifti_out)  # .dscalar.nii file
data = img.get_fdata()  # shape should be (1, 400) 
# Squeeze to remove singleton dimensions
parcel_values = np.squeeze(data)

# Show how many parcels were activated
print(f"Number of parcels with positive activation: {np.sum(parcel_values > 0)}/400")

# Extract names of parcels
def get_parcel_name(i, df_labels):
    """Extract Schaefer parcel name"""
    # i is zero-based index (0..399)
    # df_labels['Index'] is 1-based parcel ID
    parcel_id = i + 1
    # lookup the row where Index == parcel_id
    row = df_labels[df_labels['Index'] == parcel_id]
    if not row.empty:
        return row.iloc[0]['Name']
    else:
        print(f"Warning: Missing parcel {parcel_id}")
        return f"Parcel {parcel_id} (missing)"


# Extract unique network names for color coding
def extract_network(name):
    """Extract network name from Schaefer parcel name"""
    try:
        if '17Networks' in name:
            parts = name.split('_')
            if len(parts) >= 3:
                return parts[2]  # Network name
        return 'Unknown'
    except Exception as e:
        print(f"Warning: Failed to parse network from {name}: {e}")
        return 'Unknown'

# Get network names for all parcels
parcel_names = [get_parcel_name(i, df_labels) for i in range(400)]
networks = [extract_network(name) for name in parcel_names]
unique_networks = list(set(networks))

print(f"Found networks: {unique_networks}")

# Create color map for networks
colors = plt.cm.tab20(np.linspace(0, 1, len(unique_networks)))
network_color_map = dict(zip(unique_networks, colors))
bar_colors = [network_color_map[net] for net in networks]

# Improved visualization
fig, ax = plt.subplots(figsize=(14, 8))
bars = ax.bar(range(400), parcel_values, color=bar_colors, alpha=0.7)
ax.set_title('Face > Others (Colored by Network)')
ax.set_xlabel('Parcel Index')
ax.set_ylabel('Mean Contrast Value')
ax.axhline(y=0, color='black', linestyle='--', alpha=0.5)

# Add legend for networks 
handles = [plt.Rectangle((0,0),1,1, color=network_color_map[net]) for net in unique_networks]
ax.legend(handles, unique_networks, bbox_to_anchor=(1.05, 1), loc='upper left')

plt.tight_layout()
plt.show()

# Simplify labels for better readability
def simplify_label(label):
    """Simplify Schaefer labels for display"""
    try:
        # Remove '17Networks_' prefix and make more readable
        parts = label.replace('17Networks_', '').split('_')
        hemisphere = parts[0]  # LH or RH
        network = parts[1]     # Network name
        region = '_'.join(parts[2:])  # Region description
        return f"{hemisphere} {network}\n{region}"
    except:
        return label

top_n = 10
sorted_indices = np.argsort(parcel_values)[::-1]
top_indices = sorted_indices[:top_n]
top_labels = [get_parcel_name(i, df_labels) for i in top_indices]
top_values = parcel_values[top_indices]

# Extract networks for color coding
top_networks = [extract_network(label) for label in top_labels]
unique_top_networks = list(set(top_networks))
network_color_map = dict(zip(unique_top_networks, plt.cm.tab20(np.linspace(0, 1, len(unique_top_networks)))))
bar_colors = [network_color_map[net] for net in top_networks]

# Get simplified labels
simple_labels = [simplify_label(label) for label in top_labels]

# Visualize top 10 face-selective parcels
plt.figure(figsize=(14, 8))
bars = plt.bar(range(top_n), top_values, color=bar_colors, alpha=0.8, edgecolor='black', linewidth=0.5)

# Add value labels on top of bars
for i, (bar, value) in enumerate(zip(bars, top_values)):
    plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, 
             f'{value:.3f}', ha='center', va='bottom', fontsize=10)

plt.xticks(range(top_n), simple_labels, rotation=45, ha='right', fontsize=10)
plt.ylabel("Mean Contrast Value", fontsize=12)
plt.title("Top 10 Parcels: Face > Others", fontsize=14, pad=20)

# Add network legend
handles = [plt.Rectangle((0,0),1,1, color=network_color_map[net], alpha=0.8) for net in unique_top_networks]
plt.legend(handles, unique_top_networks, title='Networks', bbox_to_anchor=(1.05, 1), loc='upper left')

plt.tight_layout()
plt.show()

# Print summary table
print("\nTop 10 Face-Selective Parcels:")
print("-" * 60)
for i, (label, value, network) in enumerate(zip(top_labels, top_values, top_networks), 1):
    print(f"{i:2d}. {label:<40} | {value:6.3f} | {network}")

import os

%load_ext watermark

%watermark
%watermark --iversions

neurodesktop_version = (
    os.environ.get('JUPYTER_IMAGE', '').split(':')[-1] or
    os.environ.get('NEURODESKTOP_VERSION', 'unknown')
)

print(f"Neurodesktop version: {neurodesktop_version}")