참고 논문

fMRIflows: A Consortium of Fully Automatic Univariate and Multivariate fMRI Processing Pipelines

참고 논문 코드: https://nbviewer.org/github/miykael/fmriflows/blob/master/notebooks/02_preproc_anat.ipynb

import os
import glob
from nipype.interfaces.fsl import BET, MCFLIRT
from nipype import Node, Workflow
from nipype.interfaces.utility import Function
from matplotlib.pyplot import figure
import nibabel as nib
from nilearn.plotting import plot_stat_map

def run_bet(in_file, out_file, frac=0.5):
    from nipype.interfaces.fsl import BET
    import os

    in_file = os.path.abspath(in_file)
    out_file = os.path.abspath(out_file)

    bet = BET(in_file=in_file, out_file=out_file, mask=True, functional=True)
    bet.inputs.functional = True
    bet.inputs.frac = frac  
    result = bet.run()
    print("BET stdout:", result.runtime.stdout)
    print("BET stderr:", result.runtime.stderr)
   
    mask_file = out_file.replace('.nii.gz', '_mask.nii.gz')
    
    return out_file, mask_file

def run_mcflirt(in_file, out_file):
    from nipype.interfaces.fsl import MCFLIRT
    import os
    
    in_file = os.path.abspath(in_file)
    out_file = os.path.abspath(out_file)

    mcflirt = MCFLIRT(in_file=in_file, out_file=out_file)
    result = mcflirt.run()
    print("MCFLIRT stdout:", result.runtime.stdout)
    print("MCFLIRT stderr:", result.runtime.stderr)
    return out_file

def plot_bet_results(in_file, bet_file, bet_mask_file, plot_dir):
    import os
    import nibabel as nib  # Ensure nibabel is imported inside the function
    from nilearn.plotting import plot_stat_map
    from matplotlib.pyplot import figure
    
    img = nib.load(in_file)
    img_bet = nib.load(bet_file)
    img_mask = nib.load(bet_mask_file)
    
    fig = figure(figsize=(16, 8))
    for i, e in enumerate(['x', 'y', 'z']):
        ax = fig.add_subplot(3, 1, i + 1)
        plot_stat_map(img_bet, title=f'BET Results - {e}-axis', bg_img=img, display_mode=e, axes=ax)
    
    os.makedirs(plot_dir, exist_ok=True)
    
    base_name = os.path.basename(bet_file).replace('.nii.gz', '_bet.png')
    out_file = os.path.join(plot_dir, base_name)
    
    fig.savefig(out_file, bbox_inches='tight', facecolor='black', dpi=300, transparent=False)
    return out_file

# Define the visualization function for MCFLIRT results
def plot_mcflirt_results(in_file, mcflirt_file, plot_dir):
    import os
    import nibabel as nib  
    from nilearn.plotting import plot_stat_map
    from matplotlib.pyplot import figure
    
    img = nib.load(in_file)
    img_mcflirt = nib.load(mcflirt_file)
    
    fig = figure(figsize=(16, 8))
    for i, e in enumerate(['x', 'y', 'z']):
        ax = fig.add_subplot(3, 1, i + 1)
        plot_stat_map(img_mcflirt, title=f'MCFLIRT Results - {e}-axis', bg_img=img, display_mode=e, axes=ax)
    

    os.makedirs(plot_dir, exist_ok=True)
    
  
    base_name = os.path.basename(mcflirt_file).replace('.nii.gz', '_mcflirt.png')
    out_file = os.path.join(plot_dir, base_name)
    
    fig.savefig(out_file, bbox_inches='tight', facecolor='black', dpi=300, transparent=False)
    return out_file

base_dir = '/home/hcilab/fmri/samples' 
plot_base_directory = '/home/hcilab/fmri/plots'  

t1w_files = glob.glob(os.path.join(base_dir, '**', 'anat2', 'NIFTI', '*_T1w.nii.gz'), recursive=True)

# Check if any files are found
print(f"Found {len(t1w_files)} T1w files.")

# Process each *_T1w.nii.gz file
for t1w_file in t1w_files:
    # Ensure the t1w_file is an absolute path
    t1w_file = os.path.abspath(t1w_file)
    print(f"Processing file: {t1w_file}")
    

    subject_dir = os.path.dirname(t1w_file)
 
    sample_name = os.path.basename(os.path.dirname(os.path.dirname(subject_dir)))

    plot_directory = os.path.join(plot_base_directory, sample_name)
    
    bet_output_file = os.path.join(subject_dir, 'bet_output_file.nii.gz')
    mcflirt_output_file = os.path.join(subject_dir, 'mcflirt_output_file.nii.gz')
    
    # Create the workflow
    workflow = Workflow(name='preproc_workflow')
    workflow.base_dir = subject_dir  # Use the subject directory as the working directory

    # Create nodes for BET and MCFLIRT
    bet_node = Node(Function(input_names=['in_file', 'out_file', 'frac'],
                             output_names=['out_file', 'mask_file'],
                             function=run_bet),
                    name='bet_node')
    bet_node.inputs.in_file = t1w_file
    bet_node.inputs.out_file = bet_output_file
    bet_node.inputs.frac = 0.5  # Set this to the desired fraction

    mcflirt_node = Node(Function(input_names=['in_file', 'out_file'],
                                 output_names=['out_file'],
                                 function=run_mcflirt),
                        name='mcflirt_node')
    mcflirt_node.inputs.in_file = bet_output_file  # Use BET output as input for MCFLIRT
    mcflirt_node.inputs.out_file = mcflirt_output_file

    # Create nodes for visualization
    plot_bet_node = Node(Function(input_names=['in_file', 'bet_file', 'bet_mask_file', 'plot_dir'],
                                  output_names=['out_file'],
                                  function=plot_bet_results),
                         name='plot_bet_node')
    plot_bet_node.inputs.in_file = t1w_file
    plot_bet_node.inputs.bet_file = bet_output_file
    plot_bet_node.inputs.plot_dir = plot_directory  # Set the directory to save plots

    plot_mcflirt_node = Node(Function(input_names=['in_file', 'mcflirt_file', 'plot_dir'],
                                      output_names=['out_file'],
                                      function=plot_mcflirt_results),
                             name='plot_mcflirt_node')
    plot_mcflirt_node.inputs.in_file = t1w_file
    plot_mcflirt_node.inputs.mcflirt_file = mcflirt_output_file
    plot_mcflirt_node.inputs.plot_dir = plot_directory

    # Connect the nodes
    workflow.connect([
        (bet_node, mcflirt_node, [('out_file', 'in_file')]),
        (bet_node, plot_bet_node, [('out_file', 'bet_file'),
                                   ('mask_file', 'bet_mask_file')]),
        (mcflirt_node, plot_mcflirt_node, [('out_file', 'mcflirt_file')]),
    ])

    # Run the workflow
    workflow.run()

    # Print the output paths
    print(f"Processed {t1w_file}")
    print(f"BET results plot saved to: {os.path.join(plot_directory, os.path.basename(bet_output_file).replace('.nii.gz', '_bet.png'))}")
    print(f"Motion correction results saved to: {mcflirt_output_file}")
    print(f"MCFLIRT results plot saved to: {os.path.join(plot_directory, os.path.basename(mcflirt_output_file).replace('.nii.gz', '_mcflirt.png'))}")

참고 논문과 공통점

참고 논문 코드와 차이점