FSL预处理功能像数据

简介

FSL和AFNI以及SPM一样，是一款常用的神经影像数据处理软件。目前FSL已经出到6.0版本，我们可以通过FSL的官网来学习如何安装和使用FSL。这里是官方的安装指导页面。
要注意的是，FSL的安装是基于Python2的，并且FSL需要运行在Linux系统上。使用windows的用户可以通过虚拟机，双系统或者像我一样在移动硬盘上装一个Linux系统来使用FSL。
关于功能磁共振图像的预处理，我之前写了一篇用SPM12预处理的文章。但是由于我有用ICA-AROMA来去除运动伪迹的需求，而ICA-AROMA这个方法最好是用在FSL预处理的数据的基础上的（比如其中需要的很多文件是通过FSL预处理中自动生成的，SPM虽然也有rp开头的Realignment Parameter of time series文件，但是还是会缺少‘affmat’和/或‘warp’文件），所以我干脆又学习了一下如何用FSL来做功能磁共振图像的预处理。
ps.我后来还是用spm预处理的smooth数据做ICA-AROMA了，缺少‘affmat’和/或‘warp’文件也可以直接在标准空间做运动伪迹的矫正，具体可以看看ICA-AROMA的官方Manual。

预处理的原理

预处理的原理其实大都差不多。FSL和SPM的对比可以看看这一篇文章：fmri 分析数据 fsl & spm 两大平台比对
另外官方的预处理原理PPT文件也可以看看，同样写得非常清晰易理解。

基本步骤

前面是功能，后面是在FSL中要实现这一功能使用的函数。要查看这些函数的使用帮助，只要直接在命令行中输入函数的名字，然后回车。
其中的5~8步是要实现将功能像配准到MNI152标准空间。

1. 去除前几个时间点fslroi
2. 时间层校正 slicetimer
3. 头动校正 mcflirt
4. 去脑壳 bet
5. 功能像配准到T1加权结构像 flirt
6. T1加权结构像配准到MNI152标准空间（非线性变换）flirt
7. 得到非线性的变形场 fnirt
8. 平滑 fslmaths
9. 滤波 fslmaths
10. 每一步处理的数据都可以用fsleyes查看

具体实现

为了方便批量处理，我使用了matlab来写一个脚本文件，通过system(cmd);命令让系统调用cmd命令，从而实现一个文件批量处理预处理的多个流程。
我的数据用fsleyes的时候方向有一些问题导致normalise出来的结果很奇怪，尚不知道如何解决。这个脚本应该是没有问题的，如果有大佬路过发现问题请不吝赐教！
另外这个脚本不包括detrend和filter的部分。包含：去除时间点，Slice Timing，Realign，normalise 和 smoooth。

% FSL Preprocessing by myself
clc
clear
% FSL_path = '/usr/local/fsl';  % The Path of FSL in my PC
% spm_path = '/usr/local/MATLAB/R2019a/toolbox/spm12';
% addpath(spm_path);
cmd = 'source ~/.bashrc';
system(cmd);

%%
subjectsdir = '/home/lq/Desktop/data';  %subjects folder
subjects = {'sbj01'};
% subjects = {'sbj01','sbj02','sbj03','sbj04','sbj05','sbj07','sbj08','sbj09',...
%     'sbj10','sbj11','sbj12','sbj13','sbj14','sbj15','sbj16','sbj17','sbj18',...
%     'sbj19','sbj20','sbj21','sbj22','sbj23','sbj25','sbj26','sbj27','sbj28',...
%     'sbj29','sbj30','sbj31','sbj32','sbj33','sbj34','sbj35','sbj36','sbj38',...
%     'sbj39','sbj40','sbj41','sbj42','sbj43','sbj44','sbj45','sbj46'};
suffix = {'01','02','03','04','05','07','08','09','10',...
    '11','12','13','14','15','16','17','18','19','20',...
    '21','22','23','25','26','27','28','29','30',...
    '31','32','33','34','35','36','38','39','40',...
    '41','42','43','44','45','46'};

%%
nsubj = length(subjects);   % 被试的数量
TR = 3;
Fwhm = 6;
remove_num = 10;  % Number of timepoints removed
remain_num = 190; % Number of timepoints remained
MNI152_image_brain = '/usr/local/fsl/data/standard/MNI152_T1_2mm_brain.nii.gz';
MNI152_image = '/usr/local/fsl/data/standard/MNI152_T1_2mm.nii.gz';

for csubj = 1:nsubj
    cmd = ['mkdir ' strcat(subjectsdir,'/',subjects{csubj}) '/A'];
    system(cmd);
    cmd = ['mkdir ' strcat(subjectsdir,'/',subjects{csubj}) '/P'];
    system(cmd);
    A_func_image = strcat(subjectsdir,'/',subjects{csubj},'/','A_',suffix{csubj},'.nii.gz');
    P_func_image = strcat(subjectsdir,'/',subjects{csubj},'/','P_',suffix{csubj},'.nii.gz');
    T1w_image = strcat(subjectsdir,'/',subjects{csubj},'/','Struc_',suffix{csubj},'.nii');
    [A_cur_dir,A_func_image_name,~] = fileparts(A_func_image);
    [P_cur_dir,P_func_image_name,~] = fileparts(P_func_image);
    [Struc_cur_dir,Struc_func_image_name,~] = fileparts(T1w_image);
    A_func_image_name = A_func_image_name(1:4);
    P_func_image_name = P_func_image_name(1:4);
    
    %% 
    % Alcohol part
    % 去除时间点
    A_func_image_t = [A_cur_dir '/A/' A_func_image_name '_t.nii'];
    cmd = ['fslroi ' A_func_image ' ' A_func_image_t ' ' num2str(remove_num) ' ' num2str(remain_num)];
    system(cmd);
    
    % Slice timing
    A_func_image_ta = [A_cur_dir '/A/' A_func_image_name '_ta.nii'];
    cmd = ['slicetimer ' '-i ' A_func_image_t ' -o ' A_func_image_ta ' -r ' num2str(TR) ' --odd']; 
    % slice order=interleaved & start from odd slice (for example, from 1st slice)
    system(cmd);
    
    % Realignment
    A_func_image_tar = [A_cur_dir '/A/' A_func_image_name '_tar.nii'];
    cmd = ['mcflirt -in ' A_func_image_ta ' -refvol 0 -mats -plots ' '-o ' A_func_image_tar];
    system(cmd);
    
    % Normalise
    A_func_image_tarw = [A_cur_dir '/A/' A_func_image_name '_tarw.nii'];
    T1w_image_b = [Struc_cur_dir '/' Struc_func_image_name '_b.nii.gz'];
    A_func2T1w_mat = [A_cur_dir '/A/func2T1w_mat.mat'];
    T1w2MNI_mat = [A_cur_dir '/T1w2MNI_mat.mat'];
    T1w2MNI_warp = [A_cur_dir '/T1w_image_b'];
    
    cmd = ['bet ' T1w_image ' ' T1w_image_b ' -f 0.3 -n -m -R -B'];
%     bet <input> <output> -f 0.3 -n -m -R
    system(cmd);
    cmd = ['flirt -ref ' T1w_image_b ' -in ' A_func_image_tar ' -dof 6 -omat ' A_func2T1w_mat];
    system(cmd);
    cmd = ['flirt -ref ' MNI152_image_brain ' -in ' T1w_image_b ' -omat ' T1w2MNI_mat];
    system(cmd);
    cmd = ['fnirt --in=' T1w_image ' --aff=' T1w2MNI_mat ' --cout=' T1w2MNI_warp ' --config=T1_2_MNI152_2mm'];
    system(cmd);
    cmd = ['applywarp --ref=' MNI152_image ' --in=' A_func_image_tar ' --warp=' T1w2MNI_warp ' --premat=' A_func2T1w_mat ' --out=' A_func_image_tarw];
    system(cmd);
    
    % Smooth
    A_func_image_tarws = [A_cur_dir '/A/' A_func_image_name '_tarws.nii'];
    sigma = Fwhm/sqrt(8*log(2));
    cmd = ['fslmaths ' A_func_image_tarw ' -kernel gauss ' num2str(sigma) ' -fmean ' A_func_image_tarws];
    system(cmd);
    
    %%
    % Placebo part
    % 去除时间点
    P_func_image_t = [P_cur_dir '/P/' P_func_image_name '_t.nii'];
    cmd = ['fslroi ' P_func_image ' ' P_func_image_t ' ' num2str(remove_num) ' ' num2str(remain_num)];
    system(cmd);
    
    % Slice timing
    P_func_image_ta = [P_cur_dir '/P/' P_func_image_name '_ta.nii'];
    cmd = ['slicetimer ' '-i ' P_func_image_t ' -o ' P_func_image_ta ' -r ' num2str(TR) ' --odd']; 
    % slice order=interleaved & start from odd slice (for example, from 1st slice)
    system(cmd);
    
    % Realignment
    P_func_image_tar = [P_cur_dir '/P/' P_func_image_name '_tar.nii'];
    cmd = ['mcflirt -in ' P_func_image_ta ' -refvol 0 -mats -plots ' '-o ' P_func_image_tar];
    system(cmd);
    
    % Normalise
    P_func_image_tarw = [P_cur_dir '/P/' P_func_image_name '_tarw.nii'];
    P_func2T1w_mat = [P_cur_dir '/P/func2T1w_mat.mat'];
    
    cmd = ['flirt -ref ' T1w_image_b ' -in ' P_func_image_tar ' -dof 6 -omat ' P_func2T1w_mat];
    system(cmd);
    cmd = ['applywarp --ref=' MNI152_image ' --in=' P_func_image_tar ' --warp=' T1w2MNI_warp ' --premat=' P_func2T1w_mat ' --out=' P_func_image_tarw];
    system(cmd);
    
    % Smooth
    P_func_image_tarws = [P_cur_dir '/P/' P_func_image_name '_tarws.nii'];
    sigma = Fwhm/sqrt(8*log(2));
    cmd = ['fslmaths ' P_func_image_tarw ' -kernel gauss ' num2str(sigma) ' -fmean ' P_func_image_tarws];
    system(cmd);
end