AFNI Preprocessing & Group Analysis

import module
await module.load('afni/25.2.03') 
await module.list()

%%capture
! pip install nibabel numpy nilearn==0.13.1 scipy

import os
from glob import glob
import subprocess
import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
from nilearn import image, plotting
from ipyniivue import NiiVue
from IPython.display import display, Markdown, Image
from scipy.stats import ttest_1samp
from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor, as_completed
from tqdm.notebook import tqdm 

# 9 subjects from the Flanker Dataset
#PATTERN = "sub-0*"
# 3 subjects from the Flanker Dataset for runtime reasons
PATTERN = "sub-0[1-3]"
!datalad install https://github.com/OpenNeuroDatasets/ds000102.git
!cd ds000102 && datalad get $PATTERN

!tree -L 4 ds000102

![ -f ds000102/make_Timings.sh ] || wget -O ds000102/make_Timings.sh https://raw.githubusercontent.com/andrewjahn/AFNI_Scripts/master/make_Timings.sh

!cd ds000102 && chmod +x make_Timings.sh && bash make_Timings.sh

!cat ds000102/sub-01/func/incongruent.1D

#The following block automatically configures parallelism based on the available CPUs, 
#dividing them evenly across subjects to avoid overloading the system.
    
n_cpus = os.cpu_count() or 4 # fallback to 4 if undetectable
n_subjects = 3 # sub-01 to sub-03
max_workers = min(n_subjects, n_cpus) # don't spawn more workers than CPUs
jobs_per_subject = max(1, n_cpus // max_workers) # divide CPUs evenly across subjects

print(f"CPUs available: {n_cpus}")
print(f"Parallel subjects: {max_workers}")
print(f"Jobs per subject: {jobs_per_subject}")

# Create directory for output
os.makedirs("afni_pro_glm", exist_ok=True)

# Main subject processing function
def process_subject(subj_id, jobs=4):
    subj = f"sub-{subj_id:02d}"
    anat = f"./ds000102/{subj}/anat/{subj}_T1w.nii.gz"
    func1 = f"./ds000102/{subj}/func/{subj}_task-flanker_run-1_bold.nii.gz"
    func2 = f"./ds000102/{subj}/func/{subj}_task-flanker_run-2_bold.nii.gz"
    stim1D_c = f"./ds000102/{subj}/func/congruent.1D"
    stim1D_i = f"./ds000102/{subj}/func/incongruent.1D"
    out_dir = f"./afni_pro_glm/{subj}.results"
    script_name = f"./afni_pro_glm/proc.{subj}.tcsh"
    log_file = f"./afni_pro_glm/output.{subj}.log"

    # Skip subject if already processed
    if os.path.exists(out_dir):
        return f"{subj} skipped (already processed)"

    cmd = [
        "afni_proc.py",
        "-subj_id", subj,
        "-script", script_name,
        "-scr_overwrite",
        "-out_dir", out_dir,

        # blocks
        # NB: there is a tshift block here, but no slice timing info to apply
        "-blocks", "tshift", "align", "tlrc", "volreg", "blur", "mask", "scale", "regress",

        "-copy_anat", anat,
        "-dsets", func1, func2,

        # extra QC: radial correlation at tcat and volreg stages
        "-radial_correlate_blocks", "tcat", "volreg",

        "-tcat_remove_first_trs", "0",

        # EPI uniformization before alignment (generally helpful for human fMRI)
        "-align_unifize_epi", "local",

        # EPI-anatomical alignment options
        "-align_opts_aea",
            "-cost", "lpc+ZZ", #cost function more in line with standard starting point (and we note EPI-anatomical overlap is quite low initially)
            "-giant_move",
            "-check_flip",

        "-tlrc_base", "MNI152_2009_template_SSW.nii.gz",
        # option for anatomical-template alignment: -tlrc_opts_at -init_xform AUTO_CENTER .
        # This option is needed because the anatomical volume is very far from center to start with for anatomical-template alignment
        # + because the anatomical also has very poor overlap with the EPI dset, 
        # we might suspect that the T1w dset coordinates are not what they were when acquired.
        "-tlrc_opts_at", "-init_xform", "AUTO_CENTER",

        "-volreg_align_to", "MIN_OUTLIER",
        "-volreg_align_e2a",
        "-volreg_tlrc_warp",

        # extra QC: compute tSNR
        "-volreg_compute_tsnr", "yes",

        # set final voxel size (matching voxel grid)
        "-volreg_warp_dxyz", "3.0",

        # intersection mask of EPI and anat
        "-mask_epi_anat", "yes",

        # blur 5mm for 3x3x4 voxel data
        "-blur_size", "5.0",

        "-regress_stim_times", stim1D_c, stim1D_i,
        "-regress_stim_labels", "congruent", "incongruent",
        "-regress_basis", "GAM",
        
        # FMRI data, analyzing for voxelwise analysis
        # Motion and outlier censoring
        "-regress_censor_motion", "0.3",
        "-regress_censor_outliers", "0.05",   # extra censoring security
        "-regress_motion_per_run",

        "-regress_opts_3dD",
            "-jobs", str(jobs),
            "-gltsym", "SYM: incongruent -congruent", "-glt_label", "1", "incongruent-congruent",
            "-gltsym", "SYM: congruent -incongruent", "-glt_label", "2", "congruent-incongruent",

        "-regress_reml_exec",
        "-regress_compute_fitts",             # save fitted timeseries
        "-regress_make_ideal_sum", "sum_ideal.1D",
        "-regress_est_blur_epits",
        "-regress_est_blur_errts",
        "-regress_run_clustsim", "no",
    ]

    try:
        # Step 1: Create the AFNI processing script
        subprocess.run(cmd, check=True)

        # Step 2: Run the AFNI tcsh script, logging output
        with open(log_file, "w") as logfile:
            proc = subprocess.Popen(
                ["tcsh", "-xef", script_name],
                stdout=subprocess.PIPE,
                stderr=subprocess.STDOUT,
                text=True
            )
            for line in proc.stdout:
                print(line, end="")     # Stream to console
                logfile.write(line)     # Log to file
            proc.wait()

        return f"{subj} {'done' if proc.returncode == 0 else 'failed'}"
        
    except subprocess.CalledProcessError as e:
        return f"{subj} failed: afni_proc.py error"
    except Exception as e:
        return f"{subj} failed: unexpected error\n{e}"

# ──────────────────────────────────────────────────────────────────────────────
# Main loop to process multiple subjects in parallel
subject_ids = list(range(1, 4))  # sub-01 to sub-03
results = []
with ThreadPoolExecutor(max_workers=max_workers) as executor:
    futures = {
        executor.submit(process_subject, sid, jobs_per_subject): sid
        for sid in subject_ids
    }
    for future in tqdm(as_completed(futures), total=len(futures), desc="Subjects processed"):
        results.append(future.result())

print("\nSummary of results:")
for r in results:
    print(r)

! 3dmask_tool -input afni_pro_glm/sub-*.results/mask_group+tlrc.HEAD -prefix ./afni_pro_glm/group_mask+tlrc -inter

! 3dinfo -verb ./afni_pro_glm/sub-01.results/stats.sub-01+tlrc.BRIK

base_dir = os.path.abspath("./afni_pro_glm")

# Get all subject folders
subject_dirs = sorted(glob(os.path.join(base_dir, "sub-*.results")))

# Build paths to sub-brik 7 (e.g. Inc-Con contrast)
set_lines = []
for s in subject_dirs:
    sub_id = os.path.basename(s).replace('.results', '')
    stats_path = os.path.join(s, f"stats.{sub_id}+tlrc[7]")
    if os.path.exists(stats_path.replace('[7]', '.BRIK')) or os.path.exists(stats_path):  # extra check
        set_lines.append(f"{sub_id} {stats_path}")
    else:
        print(f"⚠️  Missing file for subject {sub_id}: {stats_path}")

# Build 3dttest++ command
cmd_list = [
    "3dttest++",
    "-prefix", os.path.join(base_dir, "group_results", "Flanker-Inc-Con_ttest"),
    "-mask", os.path.join(base_dir, "group_mask+tlrc"),
    "-setA", "Inc-Con"
]

# Add subject lines
for line in set_lines:
    cmd_list += line.split()

# Run the command
print("Running 3dttest++...")
print("\n")
print("Command:", " ".join(cmd_list))
result = subprocess.run(cmd_list, capture_output=True, text=True)

# Show output
print(result.stdout)
print(result.stderr)

# Build set list: all subjects in one set block
set_list = ["-set", "IncCon"]
for s in subject_dirs:
    sub_id = os.path.basename(s).replace(".results", "")
    coef_path = os.path.join(s, f"stats.{sub_id}_REML+tlrc[7]")
    tstat_path = os.path.join(s, f"stats.{sub_id}_REML+tlrc[8]")

    # Sanity check (optional)
    if os.path.exists(coef_path.replace("[7]", ".HEAD")) and os.path.exists(tstat_path.replace("[8]", ".HEAD")):
        set_list += [sub_id, coef_path, tstat_path]
    else:
        print(f"⚠️ Missing for {sub_id}:\n  {coef_path}\n  {tstat_path}")

# Full 3dMEMA command
cmd_list = [
    "3dMEMA",
    "-prefix", os.path.join(base_dir, "group_results", "Flanker_Inc-Con_MEMA"),
    "-mask", os.path.join(base_dir, "group_mask+tlrc"),
] + set_list

# Print and run
print("Running 3dMEMA...\n")
print("Command:", " ".join(cmd_list), "\n")
result = subprocess.run(cmd_list, capture_output=True, text=True)

print(result.stdout)
print(result.stderr)

! 3dcopy ./afni_pro_glm/group_results/Flanker-Inc-Con_ttest+tlrc ./afni_pro_glm/group_results/Flanker-Inc-Con_ttest.nii.gz
! 3dcopy ./afni_pro_glm/group_results/Flanker_Inc-Con_MEMA+tlrc ./afni_pro_glm/group_results/Flanker_Inc-Con_MEMA.nii.gz

! 3dinfo -verb ./afni_pro_glm/group_results/Flanker-Inc-Con_ttest.nii.gz

# Copy sub-brik 1 (T-statistic) for ttest++
!3dTcat -prefix ./afni_pro_glm/group_results/Flanker-Inc-Con_ttest_Tstat.nii.gz './afni_pro_glm/group_results/Flanker-Inc-Con_ttest+tlrc[1]'

# Copy sub-brik 1 (T-statistic) for MEMA
!3dTcat -prefix ./afni_pro_glm/group_results/Flanker_Inc-Con_MEMA_Tstat.nii.gz './afni_pro_glm/group_results/Flanker_Inc-Con_MEMA+tlrc[1]'

# Load T-stat maps
img_ttest = image.load_img("./afni_pro_glm/group_results/Flanker-Inc-Con_ttest_Tstat.nii.gz")
img_mema = image.load_img("./afni_pro_glm/group_results/Flanker_Inc-Con_MEMA_Tstat.nii.gz")

plotting_config = {
    "display_mode": "ortho",
    "draw_cross": False,
    "vmax": 8,
    "cmap": "cold_hot",
    "colorbar": True,
}
# transparency_range fades voxels with t < 1 to transparent
# and makes voxels with t > 3.3 (p<0.001) fully opaque
for img, title in [
    (img_ttest, "3dttest++"),
    (img_mema, "3dMEMA"),
]:
    plotting.plot_stat_map(
        img,
        title=title,
        transparency=img,
        transparency_range=[1, 3.3],
        **plotting_config,
    )

plotting.show()

! whereami -show_atlas_code

!whereami -mask_atlas_region "FS.afni.MNI2009c_asym::ctx-lh-caudalanteriorcingulate" -prefix ./afni_pro_glm/group_results/midACC_lh_mask
!whereami -mask_atlas_region "FS.afni.MNI2009c_asym::ctx-rh-caudalanteriorcingulate" -prefix ./afni_pro_glm/group_results/midACC_rh_mask

!3dmask_tool \
  -input ./afni_pro_glm/group_results/midACC_lh_mask+tlrc \
         ./afni_pro_glm/group_results/midACC_rh_mask+tlrc \
  -union \
  -prefix ./afni_pro_glm/group_results/midACC_mask

! 3dAFNItoNIFTI -prefix ./afni_pro_glm/group_results/midACC_mask.nii.gz ./afni_pro_glm/group_results/midACC_mask+tlrc

# Load AFNI mask
roi_img = image.load_img('./afni_pro_glm/group_results/midACC_mask.nii.gz')


# Plot ROI over mni background
plotting.plot_roi(roi_img,
                  title="midACC ROI Mask",
                  display_mode='ortho',
                  cmap='autumn', draw_cross=True)

! 3dresample -master ./afni_pro_glm/sub-01.results/stats.sub-01+tlrc -input ./afni_pro_glm/group_results/midACC_mask+tlrc -prefix ./afni_pro_glm/group_results/midACC_rs+tlrc

# clean up before building to ensure correct appending from scratch
for fname in [
    "./afni_pro_glm/group_results/Congruent_betas+tlrc.HEAD",
    "./afni_pro_glm/group_results/Congruent_betas+tlrc.BRIK",
    "./afni_pro_glm/group_results/Incongruent_betas+tlrc.HEAD",
    "./afni_pro_glm/group_results/Incongruent_betas+tlrc.BRIK",
]:
    if os.path.exists(fname):
        os.remove(fname)

subjects = [f"sub-{i:02d}" for i in range(1, 4)]
for subj in subjects:
    stats_path = f"./afni_pro_glm/{subj}.results/stats.{subj}+tlrc"
    
    if os.path.exists(f"{stats_path}.HEAD"):
        subprocess.run([
            "3dbucket", "-aglueto", "./afni_pro_glm/group_results/Congruent_betas+tlrc.HEAD", f"{stats_path}[1]"
        ])
        subprocess.run([
            "3dbucket", "-aglueto", "./afni_pro_glm/group_results/Incongruent_betas+tlrc.HEAD", f"{stats_path}[4]"
        ])
    else:
        print(f"Skipping {subj}: stats file not found")

# Create subject-wise contrast dataset (each sub-brick = Incongruent - Congruent)
! 3dcalc -a ./afni_pro_glm/group_results/Incongruent_betas+tlrc -b ./afni_pro_glm/group_results/Congruent_betas+tlrc -expr 'a-b' -prefix ./afni_pro_glm/group_results/incong_minus_cong+tlrc

# Extract average contrast value from ROI for each subject
! 3dmaskave -quiet -mask ./afni_pro_glm/group_results/midACC_rs+tlrc ./afni_pro_glm/group_results/incong_minus_cong+tlrc > ./afni_pro_glm/group_results/contrast_vals.txt

# Load contrast values (one per subject)
contrast_vals = np.loadtxt("./afni_pro_glm/group_results/contrast_vals.txt")

# Run a one-sample t-test against 0
t_stat, p_val = ttest_1samp(contrast_vals, 0)

# Print results
print(f"Mean contrast: {np.mean(contrast_vals):.4f}")
print(f"t = {t_stat:.4f}, p = {p_val:.4f}")

# This script creates a 5mm sphere around a coordinate
# Change the x,y,z, coordinates on the left side to select a different peak
# Radius size can be changed with the -srad option

! echo "0 20 44" | 3dUndump -orient LPI -srad 5 -master ./afni_pro_glm/group_results/Incongruent_betas+tlrc -prefix ./afni_pro_glm/group_results/ConflictROI+tlrc -xyz -

! 3dmaskave -quiet -mask ./afni_pro_glm/group_results/ConflictROI+tlrc ./afni_pro_glm/group_results/incong_minus_cong+tlrc > ./afni_pro_glm/group_results/contrast_vals_sphere.txt

# Load contrast values (one per subject)
contrast_vals_sphere = np.loadtxt("./afni_pro_glm/group_results/contrast_vals_sphere.txt")

# Run a one-sample t-test against 0
t_stat_sphere, p_val_sphere = ttest_1samp(contrast_vals_sphere, 0)

# Print results
print(f"Mean contrast: {np.mean(contrast_vals_sphere):.4f}")
print(f"t = {t_stat_sphere:.4f}, p = {p_val_sphere:.4f}")

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}")