MRI学习笔记-一阶分析

MRI数据处理需要尝试很多个不同的模型，即选择不同的时间点（onset）和持续时间（duration），所以我们会在数据文件中创建属于不同模型的文件夹，分别命名为M1、M2、M3……，下面以M1为例，后续可以更改nmodel的值。

第一步：生成该模型的onset和duration的mat文件

具体实现方法，大家可以用excel手动来弄，也可以写脚本来弄。

数据格式如图，第一列为BlockNum（session数），第二列为TrialNum（可以不用），第三列为onset_1（1号被试的onset），第四列为duration_1（1号被试的duration），第五列为condition_1（1号被试的condition），接着下面是2号被试、3号被试……

数据文件名称为onset_data_M1 （可以自己更改）

excel形式：

matlab形式：

所需要的参数都按照上述格式整理好后，需要把文件格式转换为csv，下面为matlab中将矩阵转换为csv文件。

%% onset_sort保存为csv
nmodel = 1;
T = cell2table(onset_sort(2:end,:),'VariableNames',onset_sort(1,:));
writetable(T,['onset_data_M',num2str(nmodel),'.csv']);

将excel保存为csv文件参考 怎么将excel导出为CSV格式文件-百度经验 (baidu.com)

使用matlab将csv文件读取为以第一行名称命名的矩阵

nmodel = 1;
data = readtable(['onset_data_M',num2str(nmodel),'.csv']);
names_cond = data.Properties.VariableNames;
for ai = 1:length(names_cond)
    varName = genvarname(names_cond{ai});
    eval([varName '= data.' names_cond{ai}]);
end
data = [];
save (['onset_data_M',num2str(nmodel),'.mat']);

在数据处理文件路径下新建一个文件夹，命名为 data_05_1st_level，并在该文件夹下新建文件夹，命名为 M1。将上面生成的onset_data_M1.mat复制到M1文件夹下，如下图

第二步：一阶分析前创建条件和各条件的onset和duration文件

脚本如下，可以更改的参数有

nmodel  根据自己的需要更改，1、2、3……

rootdir 数据处理最上层目录

nsess 需要处理的session个数

names 为condition名称，与条件编号顺序一致，比如我的condition为1时代表CertainGoAll，以此类推

name_sub 为被试编号

%%%一阶分析前创建条件和各条件的onset和duration文件
clc
clear
nmodel = 1;
rootdir = 'D:\motor_inhibition_fmri';%根目录
ffxdir = strcat('data_05_1st_level\M',num2str(nmodel));%结果输出目录
onsetsfile = strcat('data_05_1st_level\M',num2str(nmodel),'\onset_data_M',num2str(nmodel));%onset和duration文件位置为M1文件夹下
nsess = 12;%session个数
%% condition名称，与条件编号顺序一致
names = {'CertainGoAll','CertainStopLeft','CertainStopRight','UncertainGo',...
    'UncertainStopLeftSuccess','UncertainStopLeftFailed','UncertainStopRightSuccess','UncertainStopRightFailed',...
    'UncertainStopAllSuccess','UncertainStopAllFailed','errors'};
ncond = size(names,2);%条件个数
%% 被试编号
name_sub = 1:4;%被试编号
for ai = name_sub
    sub_no{1,ai} = num2str(ai);
end
nsub = size((sub_no),2);%计算被试数量
%% 扫描编号
for ai = 1:nsess
    session_no{1,ai} = strcat('session',num2str(ai));
end
%% 在M1文件夹中新建conditions_files
if ~exist(fullfile(rootdir,ffxdir,'conditions_files'),'dir')
    mkdir(fullfile(rootdir,ffxdir,'conditions_files'))
end
%% 读取onset.mat文件
onsetsfilename = [onsetsfile,'.mat'];
load(fullfile(rootdir,onsetsfilename));
%%
for sub = 1:nsub
    for sess = 1:nsess
        durations = cell(1,size(names,2));%生成一个1行num_condition的cell
        onsets = cell(1,size(names,2));%生成一个1行num_condition的cell
        conditems=eval(['condition_',sub_no{sub}]);%读取每个被试的条件
        onsets_sub=eval(['onset_',sub_no{sub}]);%
        duration_sub = eval(['duration_',sub_no{sub}]);%
        for i=1:ncond-1
            onsets{1,i} = onsets_sub(BlockNum == sess & conditems == i);
            durations{1,i} = duration_sub(BlockNum == sess & conditems == i);
        end
        onsets{1,ncond} = onsets_sub(BlockNum == sess & conditems == 0);
        durations{1,ncond} = duration_sub(BlockNum == sess & conditems == 0);
        %判断条件i在sess中是否不存在，不存在设为0
        isEmpty = cellfun(@isempty,onsets);
        data_nan{sub,1} = sub_no{1,sub};
        data_nan{sub,2}(sess,:) = double(isEmpty); 
        for i = 1:ncond
            if isEmpty(1,i) == 1
                onsets{1,i} = 0;
                durations{1,i} = 0;
            end
        end
        condfilename = ['conditions_',sub_no{sub},'_',session_no{sess},'.mat'];
        condfilename = fullfile(rootdir,ffxdir,'conditions_files',condfilename);
        save(condfilename,'names','onsets','durations');
    end
    fprintf('Done: %s\n',sub_no{sub})
end

第三步：建立一阶分析的batch

可更改的参数在脚本中有提示

%%建立一阶分析的batch
clc;
clear;
%% parameters
nmodel = 1;%可更改，模型编号
rootdir = 'D:\motor_inhibition_fmri';% 可更改，Directory names，最上层目录
units   = 'secs'; % 'scans' or 'secs'
imagedir = 'data_04_analysis2';%可更改，MRI数据所在文件夹
ffxdir = strcat('data_05_1st_level\M',num2str(nmodel)); % Directory to save analysis (created under subject directory)
ffxname  = strcat('ffx_M',num2str(nmodel));
rfxdir = strcat('data_06_2st_level\M',num2str(nmodel));%可更改，rfx所在文件夹
if ~exist(fullfile(rootdir,rfxdir),'dir')
    mkdir(fullfile(rootdir,rfxdir))%生成rfx的文件夹
end
conditions = {'CertainGoAll','CertainStopLeft','CertainStopRight','UncertainGo',...
    'UncertainStopLeftSuccess','UncertainStopLeftFailed','UncertainStopRightSuccess','UncertainStopRightFailed',...
    'UncertainStopAllSuccess','UncertainStopAllFailed','errors'};%可更改，conditions名称，注意顺序，后续顺序都要一致
TR      = 2; %可更改， TR in seconds slice number
tbins   = 62; %可更改， Microtime bins
tbinref = 32;% 可更改，Reference time bin for bottom-top
nsess = 12;%可更改，session的个数
name_sub = 1:4;%可更改，被试编号
for ai = name_sub
    subdir{1,ai} = num2str(ai);%被试编号
end
nsub = length(subdir);%被试数量
for ai = 1:nsess
    sessdir{1,ai}  = strcat('session',num2str(ai));
end
%% Script options
est = 1; % Include model estimation job? 1=yes
conman  = 1; % Include contrast manager job?


% Functional image file prefixes (row=subject, col=session)
functprefix = cell(nsub,nsess);
for sub = 1:nsub
    for sess = 1:nsess
        functprefix{sub,sess} = 'swa';
    end
end
%%% Only functional images to be included in the model should match functprefix in sessdir
%%% Move any others to 'dummy_scans'


% Multiple conditions files (row=subject, col=session)
multicond=cell(nsub,nsess);
for sub=1:nsub
    for sess=1:nsess
        multicond{sub,sess} = fullfile(rootdir,ffxdir,'conditions_files',['conditions_',subdir{sub},'_session',num2str(sess),'.mat']);
    end
end

% Multiple regressors (nuisance variables)
rpexp='rp_*.txt';


% Explicit mask files (row = subject)
exmask=cell(nsub,1);
for sub=1:nsub
    exmask{sub,1}=''; % Explicit mask filename
end

% Temporal / dispersion derivatives
derivs = [0 0]; % None: [0 0], time derivative: [1 0], time & dispersion: [1 1]

% High pass filter
hpf=128; % Value in seconds, or disable using 'Inf' without quotes

%% Contrasts

delcon=1; % Delete existing contrasts


% Work out contrast vector for each condition
sesscols = length(conditions)+6;
cols = sesscols*nsess+nsess;
sessind = 1:nsess;
colind = (sessind-1)*sesscols+1;
for cond = 1:length(conditions)
    thiscolind = colind+cond-1;
    eval(sprintf('%s=zeros(1,%d);',conditions{cond},cols))
    eval(sprintf('%s(thiscolind)=1;',conditions{cond}))
end

null=zeros(1,cols);

% T-contrasts
consess = [];
conpos = [];
conneg = [];
%各条件与null互相比较
% for cond = 1:length(conditions)
for cond = [1:5,7,9:length(conditions)]
    consess{1,end+1}.tcon.name = eval(['strcat(conditions{1,',num2str(cond),'} ,''-null'')']);%1:2:2*length(conditions)
    consess{1,end+1}.tcon.name = eval(['strcat(''null-'',conditions{1,',num2str(cond),'})']);%2:2:2*length(conditions)
    conpos{end+1} = conditions{1,cond};% 生成rfx需要的Positive condition
    conneg{end+1} = 'null';% 生成rfx需要的Negative condition
end
%% 条件之间比较
contrast = [1,2;1,3;1,4;2,3;4,9;5,9;7,9;9,10];%可更改，所需要比较的条件在conditions中的列数
for ncontrast = 1:length(contrast)
    a = contrast(ncontrast,1);
    b = contrast(ncontrast,2);
    consess{1,end+1}.tcon.name = strcat(conditions{1,a},'-',conditions{1,b});
    consess{1,end+1}.tcon.name = strcat(conditions{1,b},'-',conditions{1,a});
    conpos{end+1} = conditions{1,a};% 生成rfx需要的Positive condition
    conneg{end+1} = conditions{1,b};% 生成rfx需要的Negative condition
end
%生成的Positive condition和Negative condition保存到rfx的文件夹中
save(fullfile(rootdir,rfxdir,'conpos.mat'),'conpos');
save(fullfile(rootdir,rfxdir,'conneg.mat'),'conneg');
% Evaluate contrasts from names and save to text file
confile=fopen(fullfile(rootdir,ffxdir,'jobfiles','contrast_numbers.txt'),'w');
for c=1:size(consess,2)
    consess{1,c}.tcon.convec = eval(consess{1,c}.tcon.name);
    fprintf(confile,'Contrast %02d: %s\n',c,consess{1,c}.tcon.name);
end
fclose(confile);

%% Create directories, find files
ffxpath=cell(1,nsub);
scanfiles=cell(nsub,nsess);
multireg=cell(nsub,nsess);
for sub=1:nsub
    ffxpath{1,sub} = fullfile(rootdir,ffxdir,subdir{sub});
    ffxjobfilespath = fullfile(rootdir,ffxdir,'jobfiles');
    % Create ffx directory if needed

    if ~exist(ffxpath{1,sub},'dir')
        mkdir(ffxpath{1,sub});
    end

    if ~exist(ffxjobfilespath,'dir')
        mkdir(ffxjobfilespath);
    end


    for sess=1:nsess
        sesspath=fullfile(rootdir,imagedir,subdir{sub},sessdir{sess});

        % Scan filenames
        functfiles = dir(fullfile(sesspath,[functprefix{sub,sess},'*.nii']));
        if size(functfiles,1)==0
            error(['No files found using ' fullfile(sesspath,[functprefix{sub,sess} '*.nii'])])
        else
            for f=1:size(functfiles,1)
                scanfiles{sub,sess}{f,1}=fullfile(sesspath,functfiles(f,1).name);
            end
        end

        % Multiple regressor filenames
        rpfile = dir(fullfile(sesspath,rpexp));
        if size(rpfile,1)==0
            error('Multiple regressor file not found using %s', fullfile(sesspath,rpexp))
        elseif size(rpfile,1)>1
            error('More than one multiple regressor file found using %s',fullfile(sesspath,rpexp))
        else
            multireg{sub,sess} = fullfile(sesspath,rpfile.name);
        end
    end
end

%% Create jobfiles

jobfile=cell(nsub,1);
matlabbatch = cell(1,1);
for sub=1:nsub

    % Specification
    matlabbatch{1}.spm.stats.fmri_spec.dir = {ffxpath{1,sub}};
    matlabbatch{1}.spm.stats.fmri_spec.timing.units = units;
    matlabbatch{1}.spm.stats.fmri_spec.timing.RT = TR;
    matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t = tbins;
    matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t0 = tbinref;

    for sess=1:nsess
        matlabbatch{1}.spm.stats.fmri_spec.sess(sess).scans = scanfiles{sub,sess};
        matlabbatch{1}.spm.stats.fmri_spec.sess(sess).cond = struct('name', {}, 'onset', {}, 'duration', {}, 'tmod', {}, 'pmod', {});
        matlabbatch{1}.spm.stats.fmri_spec.sess(sess).multi = {multicond{sub,sess}};
        matlabbatch{1}.spm.stats.fmri_spec.sess(sess).regress = struct('name', {}, 'val', {});
        matlabbatch{1}.spm.stats.fmri_spec.sess(sess).multi_reg = {multireg{sub,sess}};
        matlabbatch{1}.spm.stats.fmri_spec.sess(sess).hpf = hpf;
    end

    matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
    matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = derivs;
    matlabbatch{1}.spm.stats.fmri_spec.volt = 1;
    matlabbatch{1}.spm.stats.fmri_spec.global = 'None';
    matlabbatch{1}.spm.stats.fmri_spec.mask = {exmask{sub,1}};
    matlabbatch{1}.spm.stats.fmri_spec.cvi = 'AR(1)';

    % Estimation
    if est==1
        matlabbatch{2}.spm.stats.fmri_est.spmmat = {fullfile(ffxpath{1,sub},'SPM.mat')};
        matlabbatch{2}.spm.stats.fmri_est.method.Classical = 1;
    end

    % Contrast manager
    if conman==1
        matlabbatch{3}.spm.stats.con.spmmat = {fullfile(ffxpath{1,sub},'SPM.mat')};
        matlabbatch{3}.spm.stats.con.consess = consess;
        matlabbatch{3}.spm.stats.con.delete = delcon;
    end

    % Save jobfile
    jobfilename = fullfile(rootdir,ffxdir,'jobfiles',[ffxname,'_',subdir{sub},'.mat']);
    jobfile{sub,1} = jobfilename;
    save(jobfile{sub,1}, 'matlabbatch')

end

% Set SPM modality back to 'fMRI' because for some reason the script
% overwrites this.
spm('defaults','FMRI');

第四步：批量运行一阶分析

可以打开spm fmri后再运行下面脚本

%%%批量运行一阶分析
clc
clear
nmodel = 1;
%每个被试batch所在的文件位置
dirroot = strcat('D:\motor_inhibition_fmri\data_05_1st_level\M',num2str(nmodel),'\jobfiles');
dirsub = strcat('\ffx_M',num2str(nmodel),'_');
name_sub = 1:4;%被试编号
for ai = name_sub
    no_sub{1,ai} = num2str(ai);
end
nsub = length(no_sub);%被试数量
for ai = 1:nsub
    fprintf('Start: %s\n',no_sub{ai});
    clear matlabbatch
    dirbatch = strcat(dirroot,dirsub,no_sub{1,ai},'.mat');
    load(dirbatch);
    spm('defaults','FMRI');
    spm_jobman('run',matlabbatch);
    fprintf('Done: %s\n',no_sub{ai});
end

至此一阶分析结束，接下来可以查看每个被试一阶分析的结果