Monday, June 16, 2014

Running SPM First-Level Analyses from the Command Line

For most people, doing first-level analyses in SPM is tedious in the extreme; for other people, these analyses are manageable through command line scripting in Matlab. And there are still other people who don't even think about first-level analyses or SPM or Matlab or any of that. This last class of people are your college classmates who majored in something else, like economics, and are living significantly less stressful lives than you, making scads of money, and regularly eating out at restaurants where they serve shrimp the size of Mike Tyson's forearm.

Assuming that you belong to the first category, however, first-level analyses are invariably a crashing bore; and, worse, completely impractical for something like beta-series analysis, where a separate regressor is needed for individual events.

The following script will do all of that for you, although it makes many assumptions about your timing files and where everything is stored. First, your directory needs to be set up with a root directory containing all of the subjects, as well as a separate timing directory containing your onsets files; and second, that the timing files are written in a four-column format of runs, regressors, onsets (in seconds), and duration. Once you have all of that, though, all it takes is simply modifications of the script to run your analysis.

As always, let me know if this actually helps, or whether there are sections of code that can be cleaned up. I will be making a modification tomorrow allowing for beta series analysis by converting all of the timing events into individual regressors, which should help things out.

Link to script: http://mypage.iu.edu/~ajahn/docs/Run_SPM_FirstLevel.m

% INPUT:
%  rootDir        - Path to where subjects are located
%  subjects       - Vector containing list of subjects (can concatenate with brackets)
%  spmDir         - Path to folder where you want to move the multiple conditions .mat files and where to the output SPM file (directory is created if it doesn't
%  exist)
%  timingDir      - Path to timing directory, relative to rootDir
%  timingSuffix   - Suffix appended to timing files
%  matPrefix      - Will be prefixed to output .mat files.
%  dataDir        - Directory storing functional runs.
%
% OUTPUT:
%  One multiple condition .mat file per run. Also specifies the GLM and runs
%  beta estimation.
%       
%
%       Note that this script assumes that your timing files are organized
%       with columns in the following order: Run, ConditionName, Onset (in
%       seconds), Duration (in seconds)
%
%       Ideally, this timing file should be written out from your
%       presentation software, e.g., E-Prime
%
%       Also note how the directory tree is structured in the following example. You may have to
%       change this script to reflect your directory structure.
%       
%       Assume that:
%         rootdir = '/server/studyDirectory/'
%         spmDir = '/modelOutput/'
%         timingDir = '/timings/onsets/'
%         subjects = [101]
%         timingSuffix = '_timing.txt'
%       1) Sample path to SPM directory:
%       '/server/studyDirectory/101/modelOutput/'
%       2) Sample path to timing file in the timing directory:
%       '/server/studyDirectory/timings/onsets/101_timing.txt'
%
% Andrew Jahn, Indiana University, June 2014
% andysbrainblog.blogspot.com

clear all

%%%----Things you will want to change for your study----%%%

rootDir = '/data/hammer/space5/ImagineMove/fmri/';
subjects = [214 215];
spmDir = '/RESULTS/testMultis/';
timingDir = 'fMRI_behav/matlab_input/';
timingSuffix = '_spm_mixed.txt';
matPrefix = 'testMulti';
dataDir = '/RESULTS/data/';

%Change these parameters to reflect study-specific information about number
%of scans and discarded acquisitions (disacqs) per run
funcs = {'swuadf0006.nii', 'swuadf0007.nii', 'swuadf0008.nii'}; %You can add more functional runs; just remember to separate them with a comma
numScans = [240 240 240]; %This if the number of scans that you acquire per run
disacqs = 2; %This is the number of scans you later discard during preprocessing
numScans = numScans-disacqs;
TR = 2; %Repetition time, in seconds

%Estimating the GLM can take some time, particularly if you have a lot of betas. If you just want to specify your
%design matrix so that you can assess it for singularities, turn this to 0.
%If you wish to do it later, estimating the GLM through the GUI is very
%quick.
ESTIMATE_GLM = 1;


%%%-----------------------------------------------------%%%



%For each subject, create timing files and jobs structure
for subject = subjects
    
    %See whether output directory exists; if it doesn't, create it
    outputDir = [rootDir num2str(subject) spmDir];
    
    if ~exist(outputDir)
        mkdir(outputDir)
    end
     
    %---Navigate to timing directory and create .mat files for each run---%
    cd([rootDir timingDir]);
    
    fid = fopen([num2str(subject) timingSuffix], 'rt');
    T = textscan(fid, '%f %s %f %f', 'HeaderLines', 1); %Columns should be 1)Run, 2)Regressor Name, 3) Onset Time (in seconds, relative to start of each run), and 4)Duration, in seconds
    fclose(fid);
    
    clear runs nameList names onsets durations sizeOnsets %Remove these variables if left over from previous analysis
    
    runs = unique(T{1});

    %Begin creating jobs structure
    jobs{1}.stats{1}.fmri_spec.dir = cellstr(outputDir);
    jobs{1}.stats{1}.fmri_spec.timing.units = 'secs';
    jobs{1}.stats{1}.fmri_spec.timing.RT = TR;
    jobs{1}.stats{1}.fmri_spec.timing.fmri_t = 16;
    jobs{1}.stats{1}.fmri_spec.timing.fmri_t0 = 1;

    %Create multiple conditions .mat file for each run
    for runIdx = 1:size(runs, 1)
            nameList = unique(T{2});
            names = nameList';
            onsets = cell(1, size(nameList,1));
            durations = cell(1, size(nameList,1));
            sizeOnsets = size(T{3}, 1);
        for nameIdx = 1:size(nameList,1)
            for idx = 1:sizeOnsets
                if isequal(T{2}{idx}, nameList{nameIdx}) && T{1}(idx) == runIdx
                    onsets{nameIdx} = double([onsets{nameIdx} T{3}(idx)]);
                    durations{nameIdx} = double([durations{nameIdx} T{4}(idx)]);
                end
            end
            onsets{nameIdx} = (onsets{nameIdx} - (TR*disacqs)); %Adjust timing for discarded acquisitions
        end

        save ([matPrefix '_' num2str(subject) '_' num2str(runIdx)], 'names', 'onsets', 'durations')
        
        %Grab frames for each run using spm_select, and fill in session
        %information within jobs structure
        files = spm_select('ExtFPList', [rootDir num2str(subject) dataDir], ['^' funcs{runIdx}], 1:numScans);

        jobs{1}.stats{1}.fmri_spec.sess(runIdx).scans = cellstr(files);
        jobs{1}.stats{1}.fmri_spec.sess(runIdx).cond = struct('name', {}, 'onset', {}, 'duration', {}, 'tmod', {}, 'pmod', {});
        jobs{1}.stats{1}.fmri_spec.sess(runIdx).multi = cellstr([outputDir matPrefix '_' num2str(subject) '_' num2str(runIdx) '.mat']);
        jobs{1}.stats{1}.fmri_spec.sess(runIdx).regress = struct('name', {}, 'val', {});
        jobs{1}.stats{1}.fmri_spec.sess(runIdx).multi_reg = {''};
        jobs{1}.stats{1}.fmri_spec.sess(runIdx).hpf = 128;
        
                
    end
    
    movefile([matPrefix '_' num2str(subject) '_*.mat'], outputDir)
    
    %Fill in the rest of the jobs fields
    jobs{1}.stats{1}.fmri_spec.fact = struct('name', {}, 'levels', {});
    jobs{1}.stats{1}.fmri_spec.bases.hrf = struct('derivs', [0 0]);
    jobs{1}.stats{1}.fmri_spec.volt = 1;
    jobs{1}.stats{1}.fmri_spec.global = 'None';
    jobs{1}.stats{1}.fmri_spec.mask = {''};
    jobs{1}.stats{1}.fmri_spec.cvi = 'AR(1)';
    
    %Navigate to output directory, specify and estimate GLM
    cd(outputDir);
    spm_jobman('run', jobs)
    
    if ESTIMATE_GLM == 1
        load SPM;
        spm_spm(SPM);
    end
    
end

%Uncomment the following line if you want to debug
%keyboard



27 comments:

  1. Hi Andy,

    As a newbie to fMRI analysis and MATLAB, your videos and tutorials have extraordinarily helpful. The library of videos that you've developed over time is truly an invaluable resource. Thank you!

    I have just a few questions about the first-level batch:

    1) Timing files: I set my paths and organized the timing files according to the structure described in the video/script comments (e.g., 01_timing.txt). However, when running the script, I ran into the following error:

    Error using textscan
    Invalid file identifier. Use fopen to
    generate a valid file identifier.

    Error in spm_first_level (line 85)

    2) Functional images (funcs): How should the functional images be named in the directory? And how should they be identified in the syntax? In your example, you reference smoothed images as 'swuadf0006.nii'. Does this represent a single image? If so, how should the rest of the series be listed? My study consists of a single run with 532 images (TR=2).

    3) I bypassed the above issues by entering timing information for a single subject (01_timing.txt, skipping textscan), and entering a single smoothed image for funcs. The script ran until line 144:

    >> spm_first_level
    Warning: Run spm_jobman('initcfg');
    beforehand
    > In spm_jobman at 107
    In spm_first_level at 144


    ------------------------------------------------------------------------
    Running job #1
    ------------------------------------------------------------------------
    Running 'fMRI model specification'

    SPM8: spm_fMRI_design (v4185) 16:01:01 - 18/06/2014
    ========================================================================
    Saving fMRI design : ...SPM.mat saved

    SPM8: spm_fmri_spm_ui (v4421) 16:01:02 - 18/06/2014
    ========================================================================
    Failed 'fMRI model specification'
    Cell contents reference from a non-cell array object.
    In file "F:\MRI_Study\dt_spm\spm8\spm_Ce.m" (v3527), function "spm_Ce" at line 53.
    In file "F:\MRI_Study\dt_spm\spm8\spm_fmri_spm_ui.m" (v4421), function "spm_fmri_spm_ui" at line 345.
    In file "F:\MRI_Study\dt_spm\spm8\config\spm_run_fmri_spec.m" (v4185), function "spm_run_fmri_spec" at line 301.

    The following modules did not run:
    Failed: fMRI model specification

    Error using cfg_util (line 835)
    Job execution failed. The full log of this
    run can be found in MATLAB command window,
    starting with the lines (look for the line
    showing the exact #job as displayed in
    this error message)
    ------------------
    Running job #1
    ------------------


    Error in spm_jobman (line 217)


    Error in spm_first_level (line 144)


    Any insight into these issues is greatly appreciated.

    Thanks,

    Roger

    ReplyDelete
  2. Hey Roger,

    Thanks for the feedback, and I'm glad to see that you're using it! It helps me debug the script and make it more accessible to other people.

    As for the questions:

    1) Check to make sure that the script is in the correct timing directory by typing "[rootDir filesep timingDir]" and seeing if the path is where you think you should be. Also test the existence of the timing file by typing "[num2str(subject) timingSuffix]" once you have specified the subject number.

    2) Each of those functional files is a full run; in your case, you only have one. Later on in the script, spm_select will grab the individual scans (or "frames" as SPM calls them) and load them into the job, which for you should be 532 images. (Remember to specify the discarded acquisitions, if any.)

    3) I would recheck steps 1) and 2) before trying to debug this; part of the script converts all of the frames into a cell, which may be getting skipped over with your manual override. In any case, include my suggestions about the first couple of steps and see what comes up.

    I'm updating the script slightly so that these errors are looked for before execution, and I am also including some additional help text to clarify things a bit. Notice that this new script will be located on my website under "Run_SPM_FirstLevel_Betas.m", since this also includes an option to do a beta series GLM. Just set the DO_BETASERIES variable to 0 to skip it.


    Best,

    -Andy

    ReplyDelete
  3. i am doing a project and have a question regarding group comparison, i have a control group and drug, if my control group is only 6 and drug is 19. I suppose it is not possible to do a group comparison am i right? if so can you further explain why

    Thank yOU

    ReplyDelete
    Replies
    1. Hi,

      It is possible to do a group comparison, but you would have to correct for the unequal sample sizes (if doing a traditional independent samples t-test). If you are using SPM and you select a two-sample t-test, the default assumptions of unequal variance and independence between groups should be appropriate. There are no restrictions against having unequal numbers of contrast or beta images in your groups.

      However, it might also be a good idea to run one-sample t-tests for each group separately and then run an ROI analysis to extract the parameters for each subject within that ROI. You could then take those values to another statistical analysis package, such as R or SPSS, and have more sophisticated correction methods at your disposal.


      Best,

      -Andy

      Delete
  4. This comment has been removed by the author.

    ReplyDelete
  5. Dear Andy,
    thank you very much for this nice script! I am new to SPM and this has certainly helped to get me started automating my pipeline. I'm not sure how to continue with the output though (always worked with FSL until now...). I think I:
    1) have to take the mean of the beta-images which represents the same condition, within-subject
    2) From this make my con-images by subtracting the right beta-images from each other, and
    3) then I can simply enter those mean images into group-level analyses.
    Does that sound about right?
    Thanks. Kind regards, Johannes

    ReplyDelete
    Replies
    1. Yep, that sounds right! You will get the average beta estimate for each regressor, which you can then convert into contrast images through the Results GUI. (I may make an automated version of that too; stay tuned).

      But yes, the rest of the steps sound right.

      -Andy

      Delete
  6. Thanks Andy, I will stay tuned then...
    FYI, I eventually came across the contrast manager (like I said, very new to SPM) and 1) just set up the right contrasts for one session and one subject, 2) using 'replicate & scale' to expand that to the other within-subject sessions, saved the .mat file and 3) finally in a for-loop just replaced the matlabbatch{1}.spm.stats.con.spmmat{1} with each subject ID. Works perfectly.
    All the best, Johannes.

    ReplyDelete
    Replies
    1. Hey Johannes,

      Nice! We use a slightly different commandline approach in our lab, but this works just as well. I would say that if you're figuring stuff like this out on your own, you're already well ahead of the pack.

      -Andy

      Delete
  7. Hi Andy,
    Great blog and script! Do you have an example of a timing textfile? I am having some troubles getting the script to read it. Thanks!

    ReplyDelete
    Replies
    1. Hey Signy,

      Yes, it would have columns like the following: 1) Run number; 2) Regressor Name; 3) Onset (in seconds, relative to start of each run); 4) Duration (in seconds)

      1 Picture 20.45 4
      1 Sound 27.00 1
      1 ButtonPress 31.42 0
      2 Sound 10.5 1
      2 Picture 15.7 4
      2 ButtonPress 22.11 0



      And so on.


      Hope this helps!

      -Andy

      Delete
    2. Thanks! And if there are multiple instances of an onset in each run, then that should be put in as a separate line?

      Delete
    3. Yes, that is correct Signy.

      -Andy

      Delete
    4. Hi Andy,
      Thanks for all your posts, videos and scripts. I am trying to run the script and I wanted to check if the Onsets should be in numerical order for each run?
      Thanks,
      Elsa Baena

      Delete
    5. Hi Elsa,

      For this script, I don't believe they need to be in numerical order. I remember running it both ordered and unordered, and the results were the same.

      -Andy

      Delete
  8. Hi again! I am getting the following error message. Any ideas why that would be? Thanks!
    SPM12: spm_fMRI_design (v5183) 11:18:57 - 16/04/2015
    ========================================================================
    Failed 'fMRI model specification'
    Error using -
    Matrix dimensions must agree.
    In file "/Users/signy/Documents/MATLAB/spm12/spm_get_ons.m" (v4855), function "spm_get_ons" at line 87.
    In file "/Users/signy/Documents/MATLAB/spm12/spm_fMRI_design.m" (v5183), function "spm_fMRI_design" at line 221.
    In file "/Users/signy/Documents/MATLAB/spm12/spm_fmri_spm_ui.m" (v6088), function "spm_fmri_spm_ui" at line 183.
    In file "/Users/signy/Documents/MATLAB/spm12/config/spm_run_fmri_spec.m" (v6157), function "spm_run_fmri_spec" at line 360.

    The following modules did not run:
    Failed: fMRI model specification

    ReplyDelete
    Replies
    1. Hi Signy,

      Usually that message means that the number of volumes in a run and the timing files aren't matching up somehow; I would doublecheck that the timing for each run (in seconds) doesn't exceed the number of TRs.

      -Andy

      Delete
  9. Hi! First, thanks for this blog. It's very helpful for new SPM users!

    I am attempting to use spm_jobman to do 2-sample t-tests for several contrasts. I saved a batch job from the gui to create my jobs structure, but when I feed this to spm_jobman, I get a very long error that seems to begin here:

    Error in spm_jobman (line 183)
    cjob = cfg_util('initjob', mljob);

    This line in jobman seems to be looking for an additional dimension in the jobs structure, but I'm not sure what belongs here, since the jobs structure created by the gui only has one dimension. Any advice would be amazing! Thanks!

    ReplyDelete
    Replies
    1. Here is the full error:

      Error using cat
      Dimensions of matrices being concatenated are not consistent.

      Error in cfg_getfile>do_filter_exp (line 1098)
      sel = unique(cat(1,ind1{:}));

      Error in cfg_getfile (line 118)
      [unused,sts] = do_filter_exp(t1,filt);

      Error in cfg_files/subsasgn_check (line 49)
      [val1, sts1] = cfg_getfile('filter',val{1},item.filter,'.*');

      Error in cfg_files/subsasgn (line 80)
      [ok, val] = subsasgn_check(item,subs,val1);

      Error in cfg_item/initialise (line 48)
      item = subsasgn(item, subs, val);

      Error in cfg_branch/initialise>initialise_job (line 54)
      item.cfg_item.val{k} = initialise(item.cfg_item.val{k}, ...

      Error in cfg_branch/initialise (line 27)
      item = initialise_job(item, val, dflag);

      Error in cfg_choice/initialise>initialise_job (line 81)
      item.cfg_item.val{1} = initialise(item.values{k}, ...

      Error in cfg_choice/initialise (line 34)
      item = initialise_job(item, val, dflag);

      Error in cfg_branch/initialise>initialise_job (line 54)
      item.cfg_item.val{k} = initialise(item.cfg_item.val{k}, ...

      Error in cfg_branch/initialise (line 27)
      item = initialise_job(item, val, dflag);

      Error in cfg_choice/initialise>initialise_job (line 81)
      item.cfg_item.val{1} = initialise(item.values{k}, ...

      Error in cfg_choice/initialise (line 34)
      item = initialise_job(item, val, dflag);

      Error in cfg_choice/initialise>initialise_job (line 81)
      item.cfg_item.val{1} = initialise(item.values{k}, ...

      Error in cfg_choice/initialise (line 34)
      item = initialise_job(item, val, dflag);

      Error in cfg_repeat/initialise>initialise_job (line 116)
      citem{l} = initialise(item.values{k}, ...

      Error in cfg_repeat/initialise (line 31)
      item = initialise_job(item, val, dflag);

      Error in cfg_util>local_initjob/@(ucjob)initialise(cjd,ucjob,false)


      Error in cfg_util>local_initjob (line 1516)
      ucj = cellfun(@(ucjob)initialise(cjd,ucjob,false), job(ui), 'UniformOutput', false);

      Error in cfg_util (line 796)
      [jobs(cjob), mod_job_idlist] = local_initjob(jobs(cjob), job, jobdedup);

      Error in spm_jobman (line 183)
      cjob = cfg_util('initjob', mljob);

      Delete
    2. Hi Tracy,

      That's a strange error, and difficult to debug on my end. If you like, I do some consulting and debugging through Google Hangouts for a fee. Email me at andrew.jahn@yale.edu if you're interested.

      Best,

      -Andy

      Delete
  10. Hello Andrew,
    I have a couple of questions.
    1. How to enter vectors of my task conditions for a number for subjects into my 1st level SPM batch.
    2. confirm how to run a for each loop to do 1st level in several subjects.
    -------
    1. I have several csv files (different file every subject) with all vectors. I would like to feed these into my SPM batch. I think the inputs of the batch go into where it says input, maybe as an array of cells, in which each input comes in the order in which it appears in the batch. To insert the values in the script is easy I think could use the function csvread (if I remember the name), and it will load in the data in matlab. Then can simply type the column of the csv that is my vector (for example :,1 (first column), :,2 (second), etc.... However, I am not sure the way I should format the vectors. Could you confirm this?

    2. So basically I think to run a for loop structured as follows:

    subject = [the names of my subjects]

    for i=1:[subject number]

    vectors = csvread(subject(i));

    [copypaste here the script of the batch]

    (...)

    input = [array made in the right format based on vectors, I don't remember the format]

    (...)

    done


    Does this make sense?

    Many many thanks
    P. Perez

    ReplyDelete
    Replies
    1. Hi Perez,

      I think that should work; just make sure that your csv files are ordered in the following columns:

      1)Run, 2)Regressor Name, 3) Onset Time (in seconds, relative to start of each run), and 4)Duration, in seconds.

      Give it a spin and let me know how it goes!

      -Andy

      Delete
  11. Hi Andy,

    Your tutorials and this blog is really helpful, especially for beginners like me. Thanks a lot.

    I have a question, while performing GLM my result's folder doesn't contain any 'beta_000k.img' file. Any idea what might have gone wrong?

    Thanks,
    Amnah

    ReplyDelete
    Replies
    1. Hi Amnah,

      Try running "Estimate," and select your SPM.mat file that was generated from setting up your first level results.

      -Andy

      Delete
  12. Hi Andy,
    I am trying to run a first level analysis with your script, but I keep getting this error message:

    Running job #1
    ------------------------------------------------------------------------
    Running 'fMRI model specification'
    Failed 'fMRI model specification'
    Error using read_hdr (line 32)
    Error reading header file "".
    In file "/Users/sheldonlab/Documents/MATLAB/spm12/@nifti/private/read_hdr.m" (v4967), function "read_hdr" at line 32.
    In file "/Users/sheldonlab/Documents/MATLAB/spm12/@nifti/nifti.m" (v4986), function "nifti" at line 26.
    In file "/Users/sheldonlab/Documents/MATLAB/spm12/spm_select.m" (v6530), function "spm_select_get_nbframes" at line 268.
    In file "/Users/sheldonlab/Documents/MATLAB/spm12/spm_select.m" (v6530), function "spm_select_expand" at line 254.
    In file "/Users/sheldonlab/Documents/MATLAB/spm12/spm_select.m" (v6530), function "spm_select" at line 130.
    In file "/Users/sheldonlab/Documents/MATLAB/spm12/config/spm_run_fmri_spec.m" (v6562), function "spm_run_fmri_spec" at line 125.

    The following modules did not run:
    Failed: fMRI model specification

    Error using MATLABbatch system
    Job execution failed. The full log of this run can be found in MATLAB command window, starting with the lines (look for the line showing the
    exact #job as displayed in this error message)

    I'm really new to Matlab, but what I get out of this is that something seems to be wrong with the read_hdr.m script, which is subsequently causing the other errors? At line 32 in read_hdr.m it says the following:

    [hdr,be] = read_hdr_raw(hname);
    if isempty(hdr)
    error('Error reading header file "%s".',hname);
    end

    The files I'm trying to process are in the Nifti format (.nii).
    I'm a little stumped what to do, doesn't that script come with the spm installation? I didn't change anything about it ... My googling has been very much unsuccessful, I'd appreciate if you could maybe clear up my confusion about what is going wrong.
    Thank you!

    Jenny

    ReplyDelete
    Replies
    1. Hi Jenny,

      That script isn't very flexible, and I'm planning to update it soon. In the meantime, I would experiment with SPM's Batch system. If you can set up a script for a single subject that is identical for all of your other subjects, then you can save the batch file and loop it over all of your subjects.


      Best,

      -Andy

      Delete
  13. Hi Andy,

    Your videos have been very helpful as I learn SPM. Thank you for taking the time to create those.

    I have a question about my timing/onset files. They are .mat files and I have a separate label for names (1 row, 12 columns), onsets (20 rows, 12 columns) , and durations (20 rows, 12 columns). I keep getting the following error when I try to run 1st level GLMS: "multiple conditions MAT-file is invalid: File must contain names, onsets, and durations cell arrays of equal length."

    I think this is because the names tab is only one row while the others are 20, but I am unsure and cannot fix it.

    Any guesses, thoughts? Thank you for your time.

    Best,
    Raquel

    ReplyDelete